AD-A172  391 


ADAPTIVE  CONTROL  WITH  FINITE  TIRE  PERSISTENCY  OF 
EXCITATIONCU)  NAVAL  POSTGRADUATE  SCHOOL  HONTEREV  CA 
I  ONUK  JUN  86 


1/2 


UNCLASSIFIED 


F/G  9/3 


■sue-  *’  VA*> 


-  - . ^ . y^yyy-  -> vvm.mi  f  r  r  r  r  ^rrf.r.Y.f.fcj 


L 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


Approved  for  public  release;  distribution  is  unlimited 


,  •>  v  v  v  v  v  vv  7  v  v  v’.  v  v  s  *.  ■  ;.  .i 


UNCLASSIFIED 

SECURITY  CLASSI^ICAl'lhfcj  OF  THIS  PAGE 


M-AI21  111 


REPORT  DOCUMENTATION  PAGE 


la  REPORT  SECURITY  CLASSIFICATION 
UNCLASSIFIED 

lb.  RESTRICTIVE  MARKINGS 

2a  SECURITY  CLASSIFICATION  AUTHORITY 

3  DISTRIBUTION/ AVAILABILITY  OF  REPORT 

Approved  for  public  release; 

Zb  DECLASSIFICATION /DOWNGRADING  SCHEDULE 

distribution  is  unlimited 

4  PERFORMING  ORGANIZATION  REPORT  NUMB£R(S) 

> 

S  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 

6a  NAME  OF  PERFORMING  ORGANIZATION 

6b  OFFICE  SYMBOL 
(If  applicable) 

7a  NAME  OF  MONITORING  ORGANIZATION 

Naval  Postgraduate  School 

Code  62 

Naval  Postgraduate  School 

6c  ADDRESS  (City,  State,  and  /IP  Code) 

7b  ADDRESS  (City,  State,  and  SIP  Code) 

Monterey,  California 


Sa  NAME  OF  FUNDING  /SPONSORING 
ORGANIZATION 


8c  ADDRESS  (C/fy,  State,  and  ZIP  Code) 


11  TITLE  ( Include  Security  Classification) 


93943-50.00 


Monterey,  California  93943-5000 


8b  OFFICE  SYMBOL 
(If  applicable) 


9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 


10  SOURCE  OF  FUNDING  NUMBERS 


PROGRAM 

PROJECT 

TASK 

WORK  UNIT 

ELEMENT  NO 

NO 

NO 

ACCESSION  NO 

ADAPTIVE  CONTROL  WITH  FINITE  TIME  PERSISTENCY  OF  EXCITATION 


1Z  PERSONAL  AUTHOR(S) 

laliWIJMWgn— — 

i3b  time  covered 

14  OATE  OF  REPORT  (Year,  Month.  Day) 

is  PAGE  COUNT 

FROM  TO 

1986,  June 

129 

■6  SUPPLEMENTARY  notation 


’  7 

COSATI  COOES  | 

18  SUBJECT  TERMS  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

F'ElD 

GROUP 

SUB-GROUP 

Adaptive  Control 

Persistency  of  Excitation 

'9  ABSTRACT  ( Continue  on  reverie  if  necessary  end  identify  by  block  number) 


An  indirect  adaptive  control  algorithm  has  been  studied  to  control 
physical  systems  with  parameter  uncertainties.  Although  the  particular 
algorithm  investigated  is  applicable  to  a  wide  class  of  discrete  time 
linear  systems,  parameter  convergence  and  therefore  global  stability  of 
the  whole  system  is  guaranteed  only  if  the  external  excitation  contains  a 
sufficiently  large  number  of  frequency  components. 


This  research  study  has  also  investigated  the  possibility  of  stopping 
the  identification  procedure  when  the  parameter  error  becomes  sufficiently 
small,  so  the  controller  automatically  turns  itself  from  adaptive  to  time 
invariant,  while  still  guaranteeing  global  stability  of  the  closed-loop 
system.  In  this  way  the  adaptive  controller  might  be  activated  only  when 
its  gains  do  not  provide  satisfactory  performances. 


.0  J  S’P'BUTlON  /  AVAILABILITY  OF  ABSTRACT 

S3  .NUASSiFiED/UNUMITED  □  SAME  AS  RPT  □  OTIC  USERS 

21  ABSTRACT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 

/.’a  SAVE  OF  RESPONSIBLE  'NOIVIOUAL 

Roberto  Cristi 

22b  TELEPHONE  (include  Area  Code) 

(408)  646-2223 

1  22 C  OFFICE  SYMBOL 

Code  62Cx 

00  FORM  1473, 84  mar 


83  APR  edition  may  be  used  until  exhausted 
All  other  editions  are  obsolete 
1 


SECURITY  CLASSIFICATION  of  *HIS  page 

UNCLASSIFIED 


V'.* 


>  >  .V 


•v\  \  N 

.SAA.V. 


'->■  '  j  w.  v 


UNCLASSIFIED 


Ml  uMli  Cl  AMi'iC  ATiom  O'  THIS  R»Ot  Sat  la«ara« 


#19  -  ABSTRACT  -  (CONTINUED) 

Computer  programs  of  the  indirect  adaptive  control 
have  been  written  for  simulation  purposes  using  both 
recursive  least  squares  and  projection  algorithms. 


'  4  ■  t  ' 


2 


SECURITY  I 


UNCLASSIFIED 


7 


MV.y-yvi; 


Approved  for  public  release;  distribution  is  unlimited. 


Adaptive  Control  with  Finite  Time 
Persistency  of  Excitation 


Irfan  Onuk 

Lieutenant  Junior  Grade,  Turkish  Navy 
B.S.,  Turkish  Naval  Academy,  1979 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


(Accession  For 

MASTER  OF  SCIENCE  IN  ELECTRICAL  ENGINEERING)  NT'"  ' s I 

|  DT*--  T  •  •'» 

!  u. 

from  the  • - - 


NAVAL  POSTGRADUATE  SCHOOL 
June  1986 


‘  :  / 


Author; 


Approved  by; 


Irf  ajv  Onuk 


Roberto  Crist i.  Thesis.  Advisor 


[.A.  Titus,  Second  Reader 


Harriett  BV  Rigas)  Chairman”  — — 
Department  of  Electrical  and  Computer  Engineering 


qm  cnrfiA, _ 

7 /  J.  N. Dyer, 

Dean  of  Science  and  Engineering 


ABSTRACT 


'  An  indirect  adaptive  control  algorithm  has  been  studied 
to  control  physical  systems  with  parameter  uncertainties. 
Although  the  particular  algorithm  investigated  is  applicable 
to  a  wide  class  of  discrete  time  linear  systems,  parameter 
convergence  and  therefore  global  stability  of  the  whole  system 
is  guaranteed  only  if  the  external  excitation  contains  a 

sufficiently  large  number  of  frequency  components. 

fk  e  •;*,  s 

This  res^aFcIr-'stttdy'' has  also  investigated  the  possibility 
of  stopping  the  identification  procedure  when  the  parameter 
error  becomes  sufficiently  small,  so  the  controller  auto¬ 
matically  turns  itself  from  adaptive  to  time  invariant,  while 
still  guaranteeing  global  stability  of  the  closed-loop  system. 
In  this  way  the  adaptive  controller  might  be  activated  only 
when  its  gains  do  not  provide  satisfactory  performances. 

Computer  programs  of  the  indirect  adaptive  control  have 
been  written  for  simulation  purposes  using  both  recursive 
least  squares  and  projection  algorithms.  -- 
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I .  INTRODUCTION 


One  of  the  most  challenging,  interesting  and  active  fields 
of  Automatic  Control  is  Adaptive  Control.  To  implement  high- 
performance  control  systems  when  the  plant  dynamics  are  poorly 
known  or  when  large  and  unpredictable  variations  occur,  the 
control  engineers  prefer  to  use  an  important  class  of  control 
systems  called  Adaptive  Control  systems.  Adaptive  Control 
comes  from  a  desire  and  need  for  improving  performance  of 
complex  engineering  systems  with  large  uncertainties.  It  is 
especially  important  in  systems  with  many  unknown  parameters 
that  are  changing  with  time.  Also,  it  can  be  defined  as  a 
special  type  of  nonlinear  feedback  control,  as  a  nonlinear, 
nonautonomous  dynamic  system. 

An  adaptive  controller  can  change  its  behavior  in  response 
to  changes  in  the  dynamics  of  the  plant  and  the  disturbances. 

The  term  adaptive  control  has  been  used  at  least  from  the 
beginning  of  the  1950's.  With  recent  advances  in  microproc¬ 
essor  technology,  it  has  become  feasible  to  implement  adaptive 
algorithms  efficiently  in  real  time  at  a  reasonable  cost. 

There  are  three  schemes  for  parameter  adaptive  control: 
gain  scheduling,  model  reference  control  and  self-tuning 
regulators.  The  starting  point  is  an  ordinary  feedback  con¬ 
trol  loop  with  a  process  and  regulator  with  adjustable 
parameters.  But  the  main  problem  is  to  find  a  convenient  way 


of  changing  the  regulator  parameters  in  response  to  changes 
in  process  and  disturbance  dynamics.  The  following  schemes 
differ  only  in  the  way  the  parameters  of  the  regulator  are 
adjusted. 

A.  GAIN  SCHEDULING 

Gain  scheduling  depends  on  finding  auxiliary  process 
variables  correlated  with  the  changes  in  plant  dynamics.  I 
this  way  it  is  possible  to  reduce  the  effects  of  parameter 
variations  by  changing  the  parameters  of  the  regulator  as 
functions  of  the  auxiliary  variables. 

The  block  diagram  of  this  scheme  is  given  in  Figure  1.1 


B.  MODEL-REFERENCE  ADAPTIVE  SYSTEMS 

In  this  type  of  adaptive  system,  a  reference  model  speci¬ 
fies  the  desired  performance,  and  tells  how  the  process  output 
should  respond  to  the  command  signal.  A  block  diagram  of 
model  reference  system  is  given  in  Figure  1.2.  As  seen,  from 
this  figure,  the  reference  model  is  part  of  the  control  system 


Figure  1.2.  Block  Diagram  of  the  Model  Reference 
Adaptive  System 


This  adaptive  system  has  two  loops :  the  lower  loop  con¬ 
tains  Regulator  and  Process.  The  upper  loop  adjusts  the 
parameters  of  the  regulator,  in  such  a  way  that  the  error 
between  the  model  output  and  the  process  output  tends  to  zero. 
In  this  type  of  control  system,  the  main  problem  is  to  deter¬ 
mine  the  adjustment  mechanism  so  that  a  stable  system  results, 
which  brings  the  tracking  error  to  zero. 

Model  reference  adaptive  systems  can  be  further  subdivided 
into  two  categories:  direct  and  indirect.  In  indirect  con¬ 
trol  the  plant  parameters  are  estimated  and  the  control  param¬ 
eters  are  adjusted  based  on  these  estimates  so  that  the  overall 
plant  transfer  function  matches  that  of  the  reference  model. 

In  direct  control  no  effort  is  made  to  identify  the  plant 
parameters  but  the  control  parameters  are  directly  adjusted  to 
minimize  the  error  between  plant  and  model  outputs. 

It  turns  out  that  the  model  reference  approach  is  applica¬ 
ble  only  to  plants  with  stable  zeroes.  In  fact  the  only  way 
the  closed  loop  transfer  function  (Plant  &  Regulator)  can 
match  the  one  of  the  model  in  its  poles  and  zeroes,  is  by  re¬ 
moving  the  plant  zeroes  by  cancellation  with  closed  loop  poles. 
This  operation  leads  to  the  presence  of  uncontrollable  or 
unobservable  modes  in  the  closed  loop  systems,  which  can  be 
accepted  only  if  they  are  stable  (i.e.,  their  effect  decays 
to  zero  with  time) ,  and  this  constitutes  the  major  limitation 
for  model  reference  adaptive  systems. 


Parameter 

estimator 


The  first  method  is  to  parameterize  the  system  directly 
in  terms  of  compensator  parameters. 

This  research  project  will  concentrate  on  the  last  method, 
that  is,  pole  placement  and  indirect  adaptive  control.  The 
evaluation  of  the  control  law  is  indirectly  determined  on 
the  basis  of  parameter  estimates.  This  method  is  also  called 
explicit,  since  the  design  is  based  on  an  explicit  estimation 
of  the  process  model..  The  parameters  of  the  controller  are 
updated  indirectly  via  estimation  of  the  process  parameters. 
Several  algorithms  to  estimate  the  parameters  of  a  linear 
model  are  available  in  the  literature.  In  this  thesic  two 
well-known  estimation  algorithms  will  be  investigated:  Recur¬ 
sive  Least  Squares  and  Projection.  They  represent  a  tradeoff 
between  complexity  of  computation  and  performances,  in  the 
sense  that  best  performances  (with  relatively  high  complexity 
are  obtained  by  using  Recursive  Least  Squares.  Also,  Block¬ 
processing  is  investigated  to  determine  the  period  of  the 
adaptation  of  the  control  parameters.  Finite  time  persistency 
of  excitation  is  also  studied  and  simulated  for  indirect 
adaptive  control . 

The  main  reason  of  choice  of  indirect  adaptive  control 
for  this  research  comes  from  the  fact  that  it  is  applicable 
to  nonminimum-phase  systems,  and  there  are  no  limitations  on 
zeroes  of  the  plant.  Therefore,  it  is  more  general  than  model 
reference  adaptive  control. 

Because  of  the  effects  of  the  zeroes  of  sampled  data 
systems  on  adaptive  control  theory,  we  start  investigating 


the  behavior  of  zeroes  of  sampled  data  systems.  This  thesis 
is  organized  as  follows:  In  Chapter  II  zeroes  of  sampled 
data  systems  are  analyzed  by  stating  results  available  in 
the  literature,  and  simulation  studies.  The  indirect  adaptive 
control  problem  is  presented  in  Chapter  III  and  computer 
simulation  studies  are  given  in  Chapter  IV. 

The  computer  programs  are  presented  in  the  Appendices  B, 

C  and  D.  Also  the  description  of  the  programs  is  given  in 
Appendix  A. 


II.  ZEROES  OF  SAMPLED  DATA  SYSTEMS 

Many  control  strategies  are  based  on  the  assumption  that 
the  zeroes  of  the  plant  are  on  a  stable  region  of  the  com¬ 
plex  plane  so  that  they  can  be  cancelled  by  precompensation 
or  closed  loop  poles.  One  example  is  the  model  reference 
adaptive  control  mentioned  in  the  Introduction.  However  it 
turns  out  that  in  sampled  data  systems  with  Zero  Order  Hold 
(the  most  popular  mean  of  Digital  to  Analog  conversion)  some 
zeroes  of  the  resulting  discrete  time  plant  are  often  in  the 
unstable  region.  In  this  chapter  we  analyze  the  position  of 
the  zeroes  of  sampled  data  systems,  in  relation  with  the  con¬ 
tinuous  time  plant  dynamics  and  sampling  frequency.  The 
structure,  the  number  of  zeroes  outside  the  unit  disc  and  the 
low  frequency  characteristic  of  the  pulse  transfer  function 
are  examined  from  the  transfer  function  for  an  important  set 
of  continuous  time  processes. 

Finally,  it  is  investigated  that  zeroes  of  the  sampled 
data  systems  are  sensitive  to  high  frequency  poles  present 
in  continuous  time  transfer  function. 

A.  TIME  DOMAIN  ANALYSIS  OF  THE  ZEROES  OF  SAMPLED  DATA 

TRANSFER  FUNCTIONS 

Poles  and  zeroes  are  important  parameters  of  linear  time- 
invariant  systems.  The  zeroes  describe  the  way  the  internal 
variables  are  coupled  to  the  inputs  and  the  outputs.  As 


described  in  [Ref.  6],  the  unstable  zeroes  limit  the  perfor¬ 
mance  of  control  systems,  since  many  design  techniques  are 
based  on  the  cancellation  of  the  process  zeroes.  This  can  be 
done  provided  the  system's  zeroes  are  stable. 

The  transformation  of  poles  in  the  continuous  time  domain 
to  the  corresponding  discrete  time  is  given  as: 

P  T 

-  e  1  (2.1) 

where  T  is  the  sampling  period.  This  transformation  maps 
the  left-half  part  of  the  s-plane  onto  the  unit  disc,  so  that 
stability  is  preserved.  But  there  is  no  simple  transformation 
for  zeroes  from  continuous  to  discrete  time  domain.  The  type 
of  hold  circuit  affects  the  position  of  the  zeroes.  Most 
digital  control  systems  use  a  zero-order  hold,  and  for  this 
type  of  hold  circuit,  the  effects  are  considered. 

In  the  following  discussion,  the  main  results  are  limit 
theorems,  which  give  the  zero  locations  for  small  and  large 
sampling  periods. 

If  the  continuous  time  transfer  function  G(S)  is  rational 
it  follows  from  Equation  (2.2)  that  H(z)  is  also  a  rational 
function 

H ( z )  =  (l-z-1)  l ( (ePiT)/(z-ePiT) )ReSpiG(s)/s  (2.2) 

i 


where : 


r\ 


P 


ReSPi 


Lim  G(s)/s(s-P.) 
s-*-Pi  1 


and  P^  are  the  poles  of  G(s)/s. 

The  function  H(z)  has  generically  n-1  zeroes.  For  particu¬ 
lar  values  of  the  sampling  period  some  zeroes  may,  however, 
go  to  infinity,  or  they  may  be  cancelled  by  poles,  i.e., 
hidden  modes.  Hidden  modes  should  not  be  considered  as  zeroes 
of  H(z).  In  fact,  there  are  in  general  no  simple  closed  form 
expressions  for  the  zeroes  of  H.  The  limiting  cases  for  small 
or  large  sampling  periods  can  be  characterized.  That  is  ex¬ 
plained  in  the  following  theorem.  The  major  steps  in  the 
proof  are  given  by  [Ref.  9] . 

Theorem  1 : 

Let  G(s)  be  a  rational  function 


(s-z-j^)  (s-z2) 
(s-p-j^)  (s-p2) 


•  (s-y 

•  (s-Pj 


(2.3) 
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and  H ( z )  the  corresponding  pulse  transfer  function.  Assume 
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that  m  <  n.  Then  as  the  sampling  period  T  -*■  0 ,  m  zeroes  of 


H(z)  go  to  1  as  exp(z^T)  and  the  remaining  n-m-1  zeroes  of 
H ( z )  go  to  the  zeroes  of  Bn_m(z)  where  Bn(z)  is  the  polynomial 
defined  as 


vTVT 
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*  .J 
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B  ( z )  =  b.n  +  b_n  +  ...  +  b 

n  I  z  2  z  n 
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The  following  results  can  be  observed  from  the  previous 


theorem: 


1.  The  limiting  zeroes  of  a  pulse  transfer  function  depend 
critically  on  the  pole  excess  of  the  corresponding 
continuous  time  system. 

2.  A  continuous  tifne  system  with  a  pole  excess  larger 
than  two  will  always  give  a  pulse  transfer  function 
with  zeroes  outside  the  unit  disc  provided  that  the 
sampling  period  is  sufficiently  short.  This  may 
happen  for  quite  reasonable  sampling  periods,  and 
sampled  data  systems  with  unstable  inverses  are  thus 
quite  common. 

An  example  is  given  to  illustrate  the  above  discussion. 

Example  2.1: 

Consider  a  system  with  transfer  function 


(s+1) 


The  corresponding  pulse  transfer  function  can  be  computed  as 


H  (z )  = 


b^z  +  b2z  +  b^ 
- 

(z-e  ) 


where : 


^  =  1  -  (1  +T  +  T2/2  )  e-T 


b2  =  (-2  +  T  +T2/2)  e  T  +  (2  +T  -T2/2)e'2T 


-  ^ 
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b3  =  (1  -T  +T2/2)e“2T  -  e-3T 

H ( z )  has  a  zero  outside  the  unit  disc  if  0  <  T  <  1.8399. 

The  result  of  Theorem  1  can  be  understood  by  studying  the 
behavior  of  the  continuous  time  transfer  function  G(s)  =  s  n 
and  its  pulse  transfer  function.  The  reason  why  Theorem  1 
holds  is  because  for  a  sampling  period  T  small  every  system 
behaves  like  G(s)  =  l/s11  ra  with  the  number  of  poles  n  and 
number  of  zeroes  m. 

The  pulse  transfer  function  corresponding  to  G(s)  =  s  n 
is  given  by 

Tn  B  ( z ) 

H  ( z)  =  - ^ - -  (2.6) 

n!  (z-l)n 

where  B^tz)  is  given  in  Equation  (2.4). 

The  polynomials  Bn  are  found  for  some  values  of  n  using 
the  program  given  in  Appendix  A. 


B1(z)  =  1 

B2  (z)  =  z  + 1 

2 

B^  (z)  =  z  +4z  +1 

B4  (z)  =  z3  +llz2  +llz  +1 

B5(z)  =  z4  +  26z3  +  66z2  +  26z  +1 

18 
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The  polynomials  Bn  have  zeroes  outside  or  on  the  unit 
circle  for  n  >_  2 .  The  unstable  zeroes  are  given  in  Table  2.1 

TABLE  2.1 

UNSTABLE  ZEROES  OF  B  (z) 

n  Unstable  Zeroes  of  Bn(z) 

2  *  -1 

3  -3.732 

4  -1,  -9.899 

5  -2.322,  -23.20 

Therefore,  it  can  be  noticed  that  there  are  continuous 
time  systems  with  stable  zeroes  such  that  the  corresponding 
pulse  transfer  function  has  unstable  zeroes. 

It  is  possible  to  give  a  complete  characterization  of  the 
zeroes  of  the  pulse  transfer  function  for  small  sampling 
periods.  A  similar  result  for  large  sampling  periods  is 
given  by  Theorem  2. 

Theorem  2 ; 

Let  G(s)  be  a  strictly  proper  rational  transfer  function 
with  G(0)  /  0  and  ReP^  <  0.  Then  all  zeroes  of  the  pulse 
transfer  function  go  to  zero  as  the  sampling  period  T  goes 


A  lower  limit  on  T  for  stable  zeroes  was  obtained  in 


[Ref.  14],  here  stated  only  for  simple  poles. 

T  >  |  £n[2A(n+l)] 


(2.7) 


where : 


6  =  max  ReP^  >  0 


(2.8) 


and 


A  =  max|Ai/G(0) 
i 


(2.9) 


If  G(0)  =  0  it  follows  from  Equation  (2.10); 


H  ( z )  =  G(0)z_1  +  (l-z-x)  l  A. 

i=l  1 


-1, 


P.T 

l 


(2.10) 


z  - 


P.T 

l 


The  corresponding  pulse  transfer  function  has  one  zero  at 
z  =  1.  The  behavior  of  the  rest  of  the  zeroes  may  be  more 
complex  as  is  shown  by  Example  2.2. 

Example  2.2: 

Let  continuous  time  transfer  function  be 


G  ( s)  = 


[  <s+ir  +1]  (s+2) 


Then  the  corresponding  pulse  transfer  function  is 


(z-1)  {  (e  i+sinT-cos  T)z+e  ‘  [1-e  ‘(sin  T+  cosT)l 


2z(z-e  )  (z  -2ze  cosT  +e  ) 


The  zeroes  are  z,  =1  and 


-2T  -T 

e  (sin  T+  cosT)  -e 

-T 

e  +  sm  T-  cos  T 


For  control  systejn  design  in  discrete  time  domain  it  is 
important  to  know  whether  the  zeroes  of  pulse  transfer  function 
are  inside  the  unit  disc  or  not.  When  all  zeroes  are  inside 
the  unit  disc  the  sampled  system  has  a  stable  inverse  and  all 
its  zeroes  can  be  cancelled.  It  is  therefore  of  interest  to 
find  sufficient  conditions  which  guarantee  that  all  zeroes  of 
a  sampled  transfer  function  are  inside  the  unit  disc.  The 
criteria  are  given  in  Theorem  3  by  [Ref.  6] . 

Theorem  3 : 


If  G(s)  is  a  strictly  proper,  rational  transfer  function 


with 


ReP .  <  0 

i 


(ii)  G  ( s)  =  0 

(iii)  -tt  <  arg  G(iw)  <  0,  for  0  <  w  <  °° 

Then  all  the  zeroes  of  the  corresponding  pulse  transfer 
function  H(z)  are  stable,  so  that  criteria  for  no  unstable 
zeroes  are  given  both  in  terms  of  G(s)  and  in  terms  of  condi¬ 
tions  on  the  Nyquist  curve  G(iw). 


B.  FREQUENCY  DOMAIN  ANALYSIS  OF  THE  ZEROES  OF  SAMPLED 

DATA  SYSTEM 

In  this  section,  the  number  of  zeroes  outside  the  unit 
disc  and  the  low  frequency  characteristic  of  the  pulse  trans¬ 
fer  function  are  analyzed  from  the  transfer  function  of  a 
set  of  continuous  time  processes. 

The  discrete  frequency  response  of  a  continuous  process 
G*(o))  can  be  derived  from  Equation  (2.11)  by  substituting 
z  =  exp(  jtoT)  : 


or  it  can  be  expressed  by  the  holding  device  transfer  function 
H(u>)  and  process  transfer  function  G (tu )  .  Then  it  becomes: 

oo 

G*(u>)  =  jp  l  H (w+kfl)  G(u>+kfl)  (2.12) 

k=-<» 


where  k  is  an  integer  and  fl  =  2tt/T.  G(oj)  is  the  frequency 
response  of  the  continuous  time  system  G(s). 

For  zero-order  hold,  the  holding  device  frequency 
response  is 

1  -e'ja)T 

H  (oj+kft)  =  ,  — - .  d-.  (2.13) 

3  (w  +  kfi ) 

Substituting  Equation  (2.13)  into  Equation  (2.12)  and  taking 
into  account  that  G(-co)  =  G(w),  where  G  is  the  conjugated 


value  of  G,  we  obtain 


G* (w )  = 


-d^t 


1  -  e 


G(w)  [1  -  f  (o)  +  e  (cj)  ] 


where : 


co  G  ( f2- co ) 
fi  -  oj  G  ( co ) 


co  v  rG(rfi+co)  G((r+1)^— oj)t  ~  ^  1(,( 

GUT  K  1  "“rn~"co  ‘  (r+1)  Q  -0)'  1  (2*16) 


[2.14; 


(2.15) 


and  r  is  an  integer. 

The  number  of  the  zeroes  outside  the  unit  disc  are  given 
by  Theorem  4  [Ref.  IQ] . 

Theorem  4 : 

If  a  continuous  process  satisfies  Equation  (2.16)  and  the 
condition 


G  ( co ) 


<  1 


(2.17) 


the  number  of  poles  of  G(z)  outside  the  unit  disc  is  zero  and 
the  phase  angle  <p(Q/2)  =  arc  G(Q/2)  is  between  0  and  -tt  then 
its  pulse  transfer  function  possesses  no  zeroes  outside  the 


unit  disc.  If  the  value  of  phase  angle  is  between  -tt  and  —  2 tt 


one  of  the  zeroes  lies  outside  the  unit  disc,  etc. 


G*(  z)  =  K* 


(z-o. 


(z-o4) 


(Z-e 


-T/T. 


( z-e 


T/T  t 


Thus,  inverse  stable  continuous  processes  frequently 
possess  inverse  unstable  pulse  transfer  functions.  For 
instance  if  in  a  transfer  function  satisfying  the  conditions 
given  by  Equations  (2.16)  and  (2.17),  t ^  >  T  and  >  T  are 

real  and  positive  for_  every  i  and  the  orders  of  its  denominator 
and  numerator  differ  by  more  than  two,  then  usually  the  pulse 
transfer  function  is  inversely  unstable. 

Also,  from  Theorem  4,  the  phase  of  the  continuous  frequency 
response  at  the  Nyquist  frequency  Q/2  is  directly  related  to 
the  number  of  unstable  discrete  time  zeroes.  High  frequency 
dynamics  (where  high  frequency  is  intended  as  above  the  Nyquist 
frequency)  influence  the  phase  at  the  Nyquist  frequency,  so 
that  it  can  cause  the  number  of  unstable  discrete  time  zeroes 
to  increase.  This  high  frequency  pole  may  come  from  the  struc¬ 
ture  of  the  measuring  devices,  actuators  or  other  hardware  parts 
of  the  physical  realization.  Usually,  we  neglect  these  terms 
in  the  transfer  function.  The  following  example  can  give  an 
idea  about  this  problem. 

Example  2.3: 

Let  the  transfer  function  of  the  plant  be 

H  ( s )  =  ~ — 

s  -1 


which  has  a  pulse  transfer  function  given  by 


r 
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•sas 


H(z)  = 


.  54308z  +  .54307 
z2  -  3.08616Z  +1 


-\'f.  ■ 

rs: 

>.  /.  / 


with  one  zero  at  z  =  -0.998. 

By  adding  the  extra  dynamics  D(s)  as 

2 


w 


D  ( s )  = 


n 


2  2 
s  +  2  ,;w„  +  w 
n  n 


the  zeroes  of  the  pulse  transfer  function  of  the  entire  sys¬ 
tem  for  c,  =  0.9  and  some  w  values  are  given  in  Table  2.2 

n  3 

and  plotted  in  Figure  2.1. 


TABLE  2.2 


'v 

ZEROES  LOCATIONS 

FOR  DIFFERENT 

w  VALUES 
n 

> 

i 

Wn 

Zeroes 

7 

V, 

7 

-3.162004 

-0.01349725 

-0.1775817 

V 

V 

V 

8 

-2.866706 

-0.009074569 

-0.1386477 

i 

9 

-2.633762 

-0.005823303 

-0.1097788 

10 

-2.447441 

-0.003552481 

-0.08844686 

% 

11 

-2.296317 

-0.002030142 

-0 .07263207 

n 

12 

-2.172057 

-0.001109288 

-0.06071956 

.■ 

V 

13 

-2.068539 

-0.0004927621 

-0.05170227 

S 

14 

-1.981240 

-0 . 0002742233 

-0.04454243 

• 

15 

-1.906793 

-0 . 00009020962 

-0.03895098 

F 

V 

V 

>.■ 

16 

-1.842641 

-0.00006002390 

-0.03435726 

17 

-1.786843 

-0.00004560289 

-0.03068041 

$ 

18 

-1.737905 

-0 . 00003143626 

-0.02745582 

v 

s' 

19 

-1.694647 

-0.00005455861 

-0.02487628 

pN 

.s 

1 

r 

20 

-1.656146 

-0.0001198984 

-0.02266212 

25 
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III.  INDIRECT  ADAPTIVE  CONTROL  FOR  DISCRETE  TIME  SYSTEM £ 


As  mentioned  in  the  Introduction,  the  main  advantage  of 
indirect  adaptive  control  systems  comes  from  the  fact  that  it 
is  applicable  to  discrete  time  systems  with  unstable  zeroes. 
It  is  assumed  that  the  system  is  as  given  in  Figure  3.1, 


Figure  3.1.  Single  Input— Single  Output  Plant  System 


p(z)  =  zn  +  P1zn  1  +  ...  +  pn  (3.2) 

r (z )  =  rxzn  1  +  r2zn~2  +  ...  +  rn  (3.3) 


The  degree  of  the  denominator  polynomial  p(z)  determines  the 
system  order  n.  Also,  if  we  assume  the  plant  to  be  causal, 
the  polynomial  r(z)  has  degree  at  most  n-1.  The  sequences 
u(t)  and  y(t)  denote  the  system  input  and  output,  respectively. 
An  alternative  way  of  describing  the  plant  is  by  its  differ¬ 
ence  equation  written  in  operator  form,  as 

p(D)y{t)  =  r (D) u (t)  (3.4) 

where : 

P (D)  =  1  +  PjD  +  ...  +  PnDn  (3.5) 

r(D)  =  r.D  +  ...  +  rnDn  (3.6) 

and  D  =  q  1  is  the  backward-shift  operator,  defined  as 
Dy  (t)  =  y (t-1)  . 

In  the  adaptive  control  problem  the  parameters  in  p(D)  and 
r(D)  are  assumed  to  be  unknown.  Only  the  system  order  n  is 
assumed  to  be  known  to  the  designer.  Hence,  two  problems 
appear  at  this  point: 

1.  Parameter  estimation;  and 

2.  Compensator  structure. 

The  general  block  diagram  of  the  indirect  adaptive  control 
problem  is  given  in  Figure  3.2.  Estimation  'f  the  parameters 
of  the  plant  is  mentioned  in  Chapter  III.C.  The  compensator 
structure  is  derived  from  the  pole  placement  problem  for 


measured.  In  general,  the  states  nay  not  be  directly  avail¬ 
able.  However,  under  certain  observability  conditions  the 
state  can  be  observed  on  the  basis  of  input-output  measure¬ 
ments.  For  this  reason  observers  are  used  to  estimate  the 
state  of  any  observable  realization.  Along  these  lines,  the 
controller  can  be  considered  as  composed  of  an  observer  and 
a  gain  matrix. 

Consider  a  state-space  representation  of  the  plant  assumed 
to  be  controllable  and  observable: 

x  ( t+1 )  =  <J>x(t)  +  Tu(t)  (3.7) 

y ( t )  =  Cx (t )  (3.8) 

with  0,  T,  C  matrices  of  dimensions  n  xn,  n  *1  and  1  xn, 
respectively. 

The  block  diagram  of  the  controlled  system  combined  with 
the  observer-controller  is  given  in  Figure  3.3.  V(t)  is  an 
external  input  which  will  be  discussed  in  Chapter  IV.  If 
we  restrict  ourselves  to  the  class  of  linear  control  inputs 
we  can  write: 


u(t)  =  -  L  x  ( t )  +v(t) 


(3.9) 


as  discussed  in  [Ref.  3],  where  L  is  a  constant  matrix  as 


L  -  [1]_  12  ...  lnJ 


•.*  */ 


vV.i 


m 
•  .■ .s  - 
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u  ( t ) 


-  Lx  (t)  +  v  (t) 


(3.12) 


To  convert  the  state-space  representation  of  the  pole- 
placement  problem  to  a  transfer  function  (using  polynomials) 
form,  taking  the  z  transform  of  Equations  (3.11)  and  (3.12) 
we  obtain 


X (z )  =  (zI-Q) _1rU (z)  +  (zI-Q) _1KY(z)  (3.13) 

U  ( z )  =  -L(zI-Q)_1rU(z)  -L  (zI-Q)  _1KY  (z)  +V(z)  (3.14) 


where  we  define  the  matrix  Q  =  $-KC.  Also  we  can  define  the 
rational  functions 


-L(zl-Q)  1r 


L  adj  (zI-Q)T 
det(zI-Q) 


k  (z) 

q(Z) 


(3.15) 


and 


-L(zl-Q)  1K 


-L  adj  (zI-Q)K 
det (zI-Q) 


...  h  (z) 
q  (z) 


(3.16) 


where  q(z),  the  observer  polynomial,  is  an  arbitrary  stable 
polynomial.  Arbitrariness  of  q(z)  is  guaranteed  by  the 
assumption  of  the  plant  being  observable,  while  h(z)  and 
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k(z),  the  controller  polynomials,  have  to  be  computed,  on  the 
basis  of  the  plant  dynamics  and  q(z). 

Thus,  the  structure  of  the  control  input  u  can  be  described 

as : 


U(z)  =  tiff  U(Z)  +  |lff  Y(Z)  +  V{Z)  (3.17) 

or  it  may  be  expressed  in  a  difference  operator  form  as 

q  (D)  u  (t)  =  K  (D)  u  ( t }  +  h(D)y(t)  +  q(D).v(t)  (3.18) 

To  make  the  notation  more  attractive,  let  us  write 
Equation  (3.18)  as 

u(t)  =  fw  u(t)  +  |w  y(t)  +  v(t)  (3-19) 

The  block  diagram  of  the  above  controlled  system  is  in  Figure 
3.4.  Hence,  the  closed-loop  system  defined  from  external 
input  v(t)  to  the  output  y(t)  has  transfer  function 

Y(t)  =  p  (D)  [q(D)-k(D?i  -r(D)h(D)  V(t)  (3.20) 

In  the  pole-placement  problem,  the  goal  is  to  determine 
the  compensator  parameters  (k,h,q)  so  that  the  closed-loop 
poles  are  as  we  desire, 


r  (D) 

*  l  r\  \ 


V  ( t ) 


y  ( t) 


(3.21) 


+ 


u(t) 


r  (D) 
p(D) 


■  .  ••  v  v 


Figure  3.4.  Transfer  Function  Form  of  the  Closed- 
Loop  System 

p*(D)  being  an  arbitrary  stable  polynomial  in  the  difference 
operator  D.  Roots  of  p*(D)  are  the  new  assigned  poles  of  the 
closed- loop  system. 

By  equating  Equations  (3.20)  and  (3.21),  we  obtain  that 
k,  h,  q  have  to  satisfy  the  polynomial  equation 

p(D)[q(D)  - k (D) ]  -  r (D) h (D)  =  q(D)p*(D)  (3.22) 

which  can  be  written  as 

k(D)p(D)  +  h(D)r(D)  =  F  ( D)  (3.23) 


where : 


F(D)  =  q (D) [p (D)  -  p* (D) ] . 
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Therefore,  the  problem  of  determining  a  suitable  compen¬ 
sator  for  pole-placement  is  equivalent  to  solving  the  poly¬ 
nomial  equation  (3.23)  which  is  known  by  the  name  of 
Diophantine  equation.  The  conditions  by  which  this  equation 
can  be  solved,  together  with  the  method  of  solution  itself, 
are  given  in  the  following  section. 

B.  DIOPHANTINE  EQUATION 

The  general  form  of  the  Diophantine  equation  in  the 
unknown  polynomials  K(z)  and  H(z)  is  given  by: 

F ( z )  =  k (z) P ( z)  +  h (z) R(z)  (3.24) 


with : 


k  ( z)  =  k  +  k.  z  +  .  .  .  +  k  zm,  k  0  (3.2  5) 

o  1  mm 

h(z)  =  h  +  h. z  +  ...  +  h  zm  (3.26) 

o  1  m 

P(z)  =  Pq  +  p^z  +  ...  +  Pnzn  (3.27) 

R(z)  =  rQ  +  r-^z  +  ...  +  r^z11  (3.28) 


and 

F ( z )  =  f  +  f.z  +  f 0z2  +  ...  +  f  zn+m  (3.29) 

o  1  2  n+m 


where : 


k.  ,  h.,  p.,  r.  and  f  arc-  constant,  not  necessarily 

all  nonzero. 

In  the  general  pole  placemen*  se*  *  .  r..: ,  the  polynomials 
p(z),  r(z),  F(z)  are  given,  »r.  z  are  unknown. 

In  the  rest  of  this  section,  we  determine  the  conditions 
by  which  the  Diophantine  equation,  can  be  solved  and  the  method 
of  solution. 

By  substituting  Equations  (  3 . 2 5 ) - ( 3 . 29 )  into  (3.24),  we 
obtain : 


c  .  x-  „  c  n+m  ,  n .  . .  . 

f  +f.  z  +  .  .  .  +  f  _  .  =  (p  +  p.  z  +  .  .  .  +  p  z  )  (k  +k.  z  +  ... 

o  1  n+m  o  1  rn  o  1 


+  k  zm)  + 
m 

(rQ  +  +  .  .  . 

+  rnzn)(ho 

h^z  +  ... 

h  z  ) 
m 

(3.30) 

and  equating  the  coefficients  of  the  same  power  of  z  yields 
the  linear  relation: 


T 

C  S 


m 


=  [f. 


"n+m 


(3.31) 


where  the  vector  C  is  defined  as: 


(3.32) 


and  the  matrix  S„  as: 


Po  Pi* 


ro  rl' 


0  p. 


0  ro- 


’  Pn-1  Pn 


,  r  ,  r 
n-1  n 


0  0 


0  0 


.p  -  p  ,  p  0 
^n-2  ^n-l 


.r  -  r  ,  r  0 
n-2  n-1  n 


0  0  .  .  .0 


0  0.  .  .0 


Equation  (3.31)  is  a  linear  algebraic  equation  in  the  unknown 
vector  C.  The  matrix  Sm  consists  of  m+1  block  rows;  each 
block  row  has  two  rows  and  can  be  obtained  by  shifting  its 
previous  block  row  to  the  right  by  one  column.  It  is  a 
2 (m+1)  x  (n+m+1)  matrix.  ' 

In  order  for  a  solution  to  exist,  the  matrix  Sm  has  to 
be  full  column  rank  [Ref.  5] .  This  can  be  satisfied  if 
2  (m+1)  .>  n+m+1,  or  m  >_  n-1. 

When  S  is  a  square  matrix  (m  =  n-1),  it  is  called  the 
m 

Sylvester  matrix  of  P  and  R,  which  has  nonzero  determinant 
provided  polynomials  P  and  R  are  mutually  coprime,  i.e.,  they 
do  not  have  common  factors  [Ref.  6]. 

We  can  summarize  the  steps  to  find  h  and  k,  mentioned 
above,  as: 

1.  Let  n  be  the  order  of  the  system. 

2.  Choose  q(z)  to  be  an  arbitrary  polynomial  of  degree  n. 


f ,  .*  -- 


3. 


Form  the  matrix  Sm  using  P(z)  and  R(z). 

4.  Set  up  the  polynomial  F(z)  as: 

F ( z )  =  q(z)[p(z)  -  p*(z)] . 

5.  Solve  for  parameters  of  h(z)  and  k(z). 

6.  Set  up  the  polynomials  h(z)  and  k(z),  or  h  (D)  and 
k(D)  using  the  relation  D  =  z"l. 

C.  PARAMETER  ESTIMATION 

As  mentioned  above,  in  indirect  adaptive  control,  we 
have  to  estimate  the  parameters  of  the  plant,  where  the 
parameters  are  the  coefficients  of  the  transfer  function 


H  ( z ) 


r  (z) 
p  (z) 


(3.34) 


We  can  compute  the  controller  parameters  h(D)  and  k(D) 
from  the  Diophantine  equation  on  the  basis  of  the  estimated 
plant  (say  r,p) . 

A  large  number  of  different  identification  methods  are 
available.  In  the  literature,  one  broad  distinction  is 
between  on-line  methods  and  off-line  methods.  In  the  off¬ 
line  case,  it  is  presumed  that  all  data  are  available  prior 
to  analysis.  Consequently,  the  data  may  be  treated  as  a 
complete  block  of  information,  with  no  strict  time  limit  on 
the  process  of  analysis.  In  contrast  to  the  off-line  case, 
the  on-line  case  deals  with  sequential  data,  which  requires 
the  paraemter  estimates  to  be  recursively  updated  within  the 
time  limit  imposed  by  the  sampling  period. 
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As  mentioned,  on-line  estimation  schemes  produce  an 
updated  parameter  estimate  within  the  time  span  between 
successive  samples. 

Also,  the  on-line  methods  are  the  only  alternative  if 
the  estimation  is  going  to  be  used  in  an  adaptive  controller 
or  if  the  process  is  time-varying.  In  this  thesis,  we  con¬ 
sider  two  classes  of  estimation  algorithms: 

1.  projection; 

2.  recursive  least-squares. 

Before  proceeding,  it  can  be  said  that  input-output 
characteristics  of  a  wide-class  of  linear  and  nonlinear 
deterministic  dynamical  systems  can  be  described  by  a  model 
expressed  in  the  following  form: 

y  (t)  =  <f>(t-l)Teo  (3.35) 


where  y(t)  denotes  the  system  output  at  time  t,  $(t-l) 
denotes  a  vector  given  as: 

lMt-l)T  =  (y(t-l)  y  ( t-2 )  ...  y(t-n)  u(t-l)  ...  u(t-n)] 

(3.36) 

and  9  denotes  the  parameters  of  the  plant,  described  as: 


9 


o 


{  ~P^ '~P2 ' ' ‘ ’ 'Pn ,rl ' r2 ' 


(3.37) 


The  following  example  illustrates  the  representation  of 


the  plant  as  given  in  Equation  (3.35). 


Example  3.1 


Suppose  the  pulse  transfer  function  is 


z  + 1 


z  +2z  +3 


It  is  also  expressed  in  the  following  form: 


H  ( z_1 ) 


-1  J  -2 
z  +  z 

-1  -2 

1  +  2  z  x  +3z  z 


-1  ^  -2 
I  +Q 


1  +2q_1  +  3q”2 


Hence,  the  difference  equation  becomes 


y ( t)  =  -2y (t-1)  -  3y ( t-2)  +  u(t-l)  +  u(t-2) 


which  can  be  expressed  as: 


=  [y (t-1)  y ( t-2)  u(t-l)  u (t-2) ] 


The  rest  of  this  section  illustrates  the  estimation 
methods.  The  projection  algorithm  and  the  recursive  least-squares 
algorithm  are  analyzed  sequentially. 


(I)  Projection  Algorithm 

By  the  Projection  Algorithm,  the  sequence  of  estimates 
0  ( t )  is  recursively  computed  as: 


=  e(t-i)  + 


_ a<t>  (t-1) _ 

c  +<Mt-l)T<Mt-l) 


-[y(t)-$(t-l)  A0  (t-1)  ]  (3.38) 


with  arbitrary  initial  estimate  0(0)  and  c  >  0;  0  <  a  <  2. 

The  detailed  derivation  steps  are  given  by  [Ref.  2]. 

This  algorithm  is  also  known  as  the  normalized  least-mean- 
squares  ( NLMS )  algorithm.  The  algorithm  results  from  the 

/v 

following  optimization  problem:  Given  0(t-l)  and  y(t0. 


determine  0(t)  so  that 


j  =  j\  |  0  (t)  -  0  ( t— 1 ) 


(3.39) 


is  minimized  subject  to 


y  ( t )  =  <J>  ( t-1 ) 1 9  (t) 


(3.40) 


The  main  properties  of  the  projection  algorithm  are  the 
following : 

1.  Estimation  of  0  at  time  t,  0 (t)  is  always  closer  to 
the  actual  value  of  0  than  the  preceding  estimated 
value  0 (t-1) . 


8(t)  -0o|  |  <  |  1 0  (t-i)  -e0 


t  >  i 


<  I  0(0)  -  0, 


:  3 . 4i ; 


2.  When  the  time  goes  to  infinity,  the  error  between  tt 
actual  output  of  the  system  and  its  predicted  value 
will  converge  to  zero. 


t-*-~  [c  +  $  (t-1)  (t-1)  ] 


where : 


(3.42) 


e ( t)  =  y(t)  -  *(t-l)  9 (t-1) 


3.  lim  |  | 0 (t)  -  0 (t-k) 


=  0  for  any  finite  k 


(3.43) 


The  adaptation  rate  converges  to  zero  as  t  -*•  °°. 


(II)  Recursive  Least  Squares  Algorithm 
According  to  Gauss  the  principle  on  which  the  least  square 
estimates  are  based  is  that  the  unknown  parameters  of  the 
process  should  be  chosen  in  such  a  way  that  the  sum  of  the 
squares  of  the  differences  between  the  actually  observed  and 
computed  values  multiplied  by  numbers  that  measure  the  degree 
of  precision  is  a  minimum. 

The  recursive  least  squares  algorithm  is  given  by  the 
following  equation  in  [Ref.  2] . 


=  0  ( t-1 )  + 


for  t  >  1  and 


.  P  ( t-2 )  $  ( t-1 )  [y ( t ) -$  (t-1)  1 0  (t-1)  ] 
1  +  4>  (t-1 )  Tp  (t-2 )  <i>  ( t-1 ) 


(3.44 


P(t-l)  =  p  (t-2)  P(t-2)^(t-l)4>(t-l)1P(t-2) 

1  +  4>  ( t-1 )  TP  ( t-2 )  $  ( t-1 ) 


(3.45) 


with  9(0)  given  and  P(-l)  is  any  positive  definite  matrix 


The  algorithm  results  from  the  minimization  of  the 


following  quadratic  cost  function: 


JN(e)  =  I  I  {y  (t)  -4>  (t-1)  T0}2  +  j(9-9(0)T  Pq-1  (0-6(0)  ) 


(3.46) 


We  can  observe  that  the  cost  function  represents  the  sum  of 
squares  of  the  output  prediction  error  e(t). 


e  (t)  =  y ( t)  -  *(t-l)10 


(3.47) 


The  second  term  of  the  right-hand  side  of  the  cost  function 
accounts  for  the  initial  parameter  estimates  weighted  by  the 
matrix  PQ.  In  this  way,  we  can  consider  PQ  (the  initial  con¬ 
dition  of  P(t)  in  Equation  (3.45))  as  the  "confidence"  on 
the  initial  conditions  0(0). 

The  least  squares  algorithm,  as  can  be  seen  from  the 
simulation  results,  has  much  faster  convergence  than  the 
projection  algorithm. 


D.  PERSISTENCY  OF  EXCITATION 

In  order  to  guarantee  global  stability  of  indirect  adap¬ 
tive  control  algorithms,  the  input  to  the  plant  has  to  be 
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persistently  exciting,  which  in  turn  implies  global  conver¬ 
gence  of  the  plant  parameters  to  the  corresponding  actual 
values  [Ref.  3] .  To  obtain  consistent  estimates  of  the  plant 
parameters,  it  is  necessary  that  the  input  signal  to  the  plant 
be  sufficiently  rich  in  frequencies,  and  excites  all  modes 
of  the  plant. 

In  a  closed  loop  set-up,  the  control  input  Is  the  sum  of 
an  external  input  signal  v(t)  and  a  feedback  signal  from 
the  adaptive  controller. 

The  feedback  signal  may  in  principle  cancel  any  excitation 
contained  in  the  external  input  signal.  This  problem  and 
the  potential  for  unbounded  growth  of  the  control  and  output 
signals  have  made  the  guarantee  of  persistent  excitation  a 
difficult  problem.  Recently,  Elliot  [Ref.  8]  gave  sufficient 
conditions  which  guarantee  persistency  of  excitation.  The 
following  theorems  summarize  these  conditions. 

Theorem  4.1: 

Let  w  be  the  number  of  parameters  estimated.  In  order  to 
guarantee  global  convergence  of  the  plant  parameters  to  their 
true  values,  the  following  conditions  have  to  be  satisfied: 

1.  The  external  input  v(t)  should  consist  of  a  sum  of  2w 
sinusoids . 

2.  The  compensator  parameters  should  be  updated  each  N 
samples,  with  N  >_  lOn. 

Therefore,  in  this  thesis  report,  block  processing  is 
used  in  the  sense  that  N  data  samples  are  taken,  and  N  iter¬ 
ations  of  the  recursive  least-squares  algorithm  are  performed 
between  control  parameters  (h,k)  updates.  To  avoid  time 
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variation  of  the  plant  dynamics  during  the  spanning  process, 
the  control  parameters  are  held  constant  during  spanning 
blocks  (of  length  N)  and  are  changed  only  between  them. 

As  we  can  see  from  Theorem  4.1,  these  two  conditions 
are  sufficient  to  guarantee  parameter  convergence  for  any 
initial  conditions  as: 

A 

1 im  0 ,  • =  0  * 

k 

A 

being  the  estimate  of  0*  at  time  t^. 

In  the  next  section,  we  will  examine  finite  time  persis¬ 
tency  of  excitation. 

E.  FINITE  TIME  PERSISTENCY  OF  EXCITATION 

It  is  clear  that  the  persistency  of  excitation  condition 
on  the  external  input  v(t)  is  the  main  limitation  about  this 
adaptive  control  algorithm.  Recently,  in  a  report  by  Crist i 
[Ref.  15],  a  possible  solution  to  this  problem  has  been  given, 
by  stopping  parameter  adaptation  when  the  performances  are 
close  to  the  desired  ones.  More  precisely,  the  model  output 
error  between  desired  output  of  the  closed-loop  system  and 
controlled  plant's  output  is  the  measure  of  how  far  the  system 
is  from  the  desired  performance  (i.e.,  pole  placement)  and 
adaptation  is  stopped  whenever  the  error  falls  below  an  arbi¬ 
trary,  preassigned  threshold.  The  configuration  of  the  system 
is  given  in  Figure  3.5. 


A  major  difficulty  with  this  algorithm  is  to  be  able  to 
guarantee  global  stability  of  the  whole  system.  Therefore 
we  have  to  prove  that  the  signals  in  the  loop  are  bounded, 
provided  the  external  input  v(t)  is  bounded.  Global  stability 
of  the  closed-loop  system  is  proved  below. 

When  time  goes  to  infinity  and  persistency  of  excitation 
is  present,  parameters  of  the  estimates  of  the  plant  poly- 
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nomials  rfc  and  P  converge  to  the  actual  values  r  and  p. 
Therefore,  the  difference  between  measured  output  and  desired 
output  tends  to  zero,  and  also  U(t)  tends  to  zero.  Thus,  by 
definition  of  limit,  there  exists  an  instant  tQ  such  that 

I u (t)  |  <  e 

for  all  t  >  t  and  for  any  given  c  >  0.  Rearranging  Equation 
(3.48)  yields 

y (t)  =  § v(t)  +  u{t)  (3-49) 

which  can  be  expressed  on  a  block  diagram  form  as  in  Figure  3.6 


Figure  3.6.  Block  Diagram  Representation  of  Equation 


It  can  be  interpreted  as  a  single  input  single  output 
linear  stable  process  with  bounded  input  v(t)  and  bounded 
disturbance  u(t).  From  the  fact  that  stable  linear  systems 
with  a  bounded  input  produce  a  bounded  output,  the  global 
stability  results  follow  easily.  So  that  with  finite  time 
persistency  excitation,  the  closed-loop  system  is  eventually 
equivalent  to  a  system  with  poles  as  desired,  bounded  zeroes 
with  a  bounded  output,  disturbance  y(t). 

We  can  generate  the  external  input  v(t)  as 

v(t)  =  ve(t)  +  vc(t)  (3.50) 

with  v_(t)  the  external  desired  command  and  v  (t)  an  added 
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persistency  excitation  signal.  When  the  adaptation  is  con¬ 
tinuing,  the  external  input  is  composed  of  v  (t)  and  v  (t) , 
but  after  stopping  the  adaptation  persistency  of  excitation 
is  not  needed,  it  will  be  identical  to  the  longer,  so  that 
vc(t)  can  decay  to  zero. 

In  this  way,  the  adaptive  algorithm  is  activated  only 
when  the  performance  index  u(t)  is  larger  than  a  minimum 
threshold.  In  the  next  chapter,  simulation  studies  will  show 
their  efficiency  via  some  examples. 


SIMULATION  STUDIES 


This  research  report  includes  three  computer  programs. 

The  program  named  CONDIS,  given  in  Appendix  B,  is  used  to 
investigate  the  behavior  of  the  zero  in  sampled  data  systems. 
This  program  computes  the  parameters  of  the  discrete  time 
transfer  function  from  the  parameters  of  its  continuous  time, 
and  the  sampling  rate.  This  program  has  been  used  to  inves¬ 
tigate  the  perturbations  of  the  zeroes  for  different  values 
of  the  sampling  interval. 

The  other  two  programs  given  in  Appendix  C  and  D  simulate 
the  indirect  adaptive  control  using  recursive  least  squares 
and  projection  algorithms.  Entering  order  and  numerator, 
denominator  parameters  of  the  plant,  observer  polynomial  and 
desired  closed-loop  characteristic  polynomial,  the  program 
simulates  the  entire  system  in  an  interactive  fashion.  Also 
it  is  capable  of  graphical  and  tabulation  results. 

The  adaptive  controller  presented  above  has  been  simulated 
for  several  different  plant  dynamics. 

Example  4.1: 

Let  the  discrete  time  transfer  function  of  the  plant  be: 


z  +  2 _ 

z2  -  2z  +0.75 


z  +  2 

(z  -1.5)  (z  -0.5) 


The  plant  has  one  unstable  zero,  z  =  -2,  and  two  poles 


p.  =  1.5  and  p0  =  0.5. 


The  polynomials  q(z)  and  p*(z)  are  chosen  to  be  stable 


polynomials  with  degree  n. 

In  particular,  let 

q  (z )  =  z2  -  0.3z  -  0.28  =  (z-0.7)(z+0.4) 


and 


.  P* (z) 


2 

z 


l.lz  +  0.3 


(z  -0.5)  (z  -0.6) 


both  having  stable  roots,  i.e.,  inside  the  unit  disc  in  the 
z-plane.  Before  going  into  simulation  study,  the  problem  is 
solved  analytically  by  comparing  their  responses.  By  writing 
the  plant  dynamics  in  shift-operator  form. 


H  (q 


-1, 


-iT.-tlgl2, 


1  - 2q_1  +0.75q'2 


we  can  derive  the  difference  equation  of  the  given  system  as 


y ( t )  =  2y ( t-1 )  -  0 . 75y (t-2)  +  u(t-l)  +  2u(t-2) 

y  ( t)  =  <t>(t-l)Te* 


where 


and 


u(t-l) 
u  ( t-2 ) 


* 

0 


-2 

0.75 

1 

2 


0*  being  the  vector  of  the  actual  parameters  of  the  plant. 
Generally,  it  can  be  partitioned  as: 


0*  = 


£  and  r  being  vectors  of  parameters  of  the  denominator  and 
numerator  of  the  plant,  respectively. 

From  Equation  (3.23),  considering  k(z)  and  h(z)  as 
controller  polynomials,  we  can  write 


k  ( z )  p  ( z)  +  h(z)r(z)  =  q(z)[p(z)  -p*(z)] 


In  this  example  we  have  chosen: 

2 
z 


3$ 


I 


p(z) 


2z  +  0.75 


ffl 


p*(z)  =  z  -l.lz+0.3 


By  the  constraints  on  the  polynomials  for  solvability  of 
the  Diophantine  equation,  we  choose  k  and  h  to  be  first  order 
polynomials  as: 


k(z)  =  k^z  +  kQ 


h (z)  =  h1z  +  hQ 


Substituting  these  polynomials  into  the  Diophantine 
equation  in  matrix  form  yields: 


0.75  -  2  1  0 


-0.126 


hQ  2  10  0 


0.117 


0  0.75  -2  1 


0.72 


2  10 


-0.9 


Since  the  plan  does  not  have  any  pole-zero  cancellation, 


the  determinant  of  the  matrix  S  is  not  equal  to  zero.  Hence, 
it  is  solvable.  Solving  the  above  matrix  equation,  we  compute 
the  control  parameters  k^  =  -.899,  kQ  =  -0.688,  h^  =  -0.3900 
and  hQ  =  0.195.  The  corresponding  polynomials  of  the  con¬ 
troller  are: 


k  (z )  =*  -  .  899z  -  0.688 

h (z)  =  0 . 39z  +  0.195 


This  corresponds  to  the  optimal  choice  of  compensator 
parameters  for  the  desired  pole  placement.  Therefore  the 
compensator  given  by  the  difference  equation 


u(t) 


k  ( D )  u  ( t)  +  h  (D)  y  (t ) 
q(D) 


+  v  ( t ) 


yields  closed-loop  transfer  function 


H  ( z ) 


Y(z) 

vTzT 


_ q(z)r  (z) _ 

p  ( z )  [q  ( z )  -  k  ( z )  ]  -  r  (z )  h  ( z ) 


which  becomes: 


H(z) 


z  +2 

z2  —  1 . 1  z  +0.3 


.  V*  •/ 


The  step  response  of  this  system  is  given  in  Figure  4.1. 

The  next  approach  to  the  problem  is  simulation  of  the  system 
using  a  computer  program.  It  is  observed  that  from  the  simu¬ 
lations,  after  the  lOn  iterations,  the  estimated  parameter 
will  be  closer  to  0*  and  0  will  tend  to  zero,  where 

0  =  ||e*  - e  J  | 

A  A 

This  means  that  p  and  r  are  converging  to  the  actual  param- 
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eters,  and  h  and  k  converge  to  the  actual  values  for  the 
desired  pole-placement.  The  desired  output  and  controlled 
system's  output  are  given  in  Figure  4.2.  Also,  parameter 
error  and  output  prediction  error  are  given  in  Figures  4.3 
and  4.4,  respectively.  After  convergence  of  the  prediction 
error  to  zero,  the  persistency  of  excitation  added  to  the 
external  command  is  turned  off.  In  this  example,  we  have 
chosen  the  step  as  external  command.  External  input  and  plant 
input  are  in  Figures  4.5  and  4.6,  respectively. 

If  we  compare  Figures  4.1  and  4.2,  analytical  and  simu¬ 
lated  system's  results  are  almost  the  same. 

In  the  next  example,  it  is  considered  that  one  of  the 
plant  parameter  is  changed  after  some  period  of  time.  It  is 
investigated  how  the  system  behavior  will  change  in  this 
condition . 

Example  4.2: 


y.vS 

-  «.*  vv, 


f.v^s 


VV 
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Let  the  plant  be  the  same  as  the  one  in  the  previous 
example.  After  two  blocks  length,  the  zero  of  the  plant  is 


0005 


PARAMETER  ERROR 


OUTPUT  PREDICTION  ERROR 
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EXTERNAL  INPUT 


INPUT  TO  THE  PLANT 


Input  vs 


perturbed  from  z  =  -2  to  z  =  -3.  To  compute  the  parameters 
associated  to  the  perturbed  plant  with  pulse  transfer  function 

H  ( z)  =  -= - -"+-- - 

z  -  2z  +  0.75 

write  its  difference  equation  as 


y(t)  =  2y (t-1)  -  0 . 75y (t-1)  +  u(t-l)  +  3u(t-2) 


and  immediately 

-2.0 

0.75 

6*  = 

1.0 

3.0 


By  using  the  same  observer  and  desired  polynomials  g(z), 

p*(z)  and  following  the  same  steps  as  in  Example  4.1,  the 

controller  parameters  become  k^  =  -0.894,  kQ  =  -0.778, 

h.  =  -0.305  and  h  =  0.152. 

1  o 

Then,  the  polynomials  k(z)  and  h(z)  are: 


k (z)  =  - .  89 4 z  -  0.778 


h (z)  =  -0.305z  +  0.152 


The  closed  loop  transfer  function  becomes: 


z  -  1.12  +  0.3 

.  The  step  response  of  this  transfer  function  is  given  in 
Figure  4.7.  Also,  the  other  graphical  results  are  in  Figures 
4.8-4.12. 

In  the  next  example,  we  investigate  how  the  disturbance 
at  the  output  will  affect  the  behavior  of  the  controlled 
system. 

Example  4.3: 

Consider  a  plant  with  pulse  transfer  function  as: 


n\£.i  —  — x - 

z  -  2z  +  0.75 

and  assume  an  output  disturbance  exists.  Let  the  disturbance 
be  sinusoidal  with  frequency  5tt/2  rad/sec.  In  this  case,  it 
is  observed  that  the  estimated  parameters  of  the  plant  don't 
converge  to  the  actual  parameters.  Corresponding  plots  are 
given  in  Figure  4.13  through  4.17.  However,  if  the  disturbance 
is  small,  the  parameters  converge  close  to  the  actual  values 
and  we  can  still  obtain  satisfactory  performances. 

Comparison  of  Different  Estimation  Techniques:  RLS  and  P.A. 

In  the  above  examples  we  used  a  recursive  least-squares 
algorithm  for  estimating  the  plant  parameters.  In  this  sec¬ 
tion,  we  compare  the  behavior  of  the  indirect  adaptive  control, 
when  the  projection  algorithm  is  used. 


The  projection  algorithm  has  the  advantage  of  requiring 
less  computations.  Approximately,  the  number  of  computations 
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grows  as  n  in  the  recursive  least-squares  algorithm.  But 


in  the  projection  algorithm,  it  grows  as  n.  As  noted  before. 


the  projection  algorithm  can  be  described  as 


=  0  ( t-1 )  + 


a$ (t-1) 


c  +4  (t-1)  4  (t-1) 


•[y  (t)  -  4  (t-1 ) 1 0  (t-1 )  ] 


with  c  >  0  and  0  <  a  <  2.  The  value  of  the  parameter  a  is 


crucial.  The  behavior  os  the  algorithm  greatly  depends  on 


its  value. 


In  the  next  example,  the  projection  algorithm  is  used 


for  estimation  procedure  on  the  previous  process. 


Example  4.4: 


Using  the  same  pulse  transfer  function  and  controller 


polynomials  as  in  Example  4.3  and  choosing  the  constants 


a  =  1  and  c  =  1  in  the  above  equation,  results  can  be  observed 


from  Figure  4.18  through  4.22. 


Actually,  using  this  algorithm,  several  systems  have 


been  simulated.  This  result  is  the  most  reasonable  one. 


Considering  the  above  results,  indirect  adaptive  control 


is  a  very  effective  control  technique  when  the  parameters  of 


the  plant  are  unknown  and  large  and  unpredictable  variations 


occur  in  the  plant  dynamics.  Also,  it  is  applicable  to 


nonminimum- phase  systems.  This  can  be  considered  as  a  big 


advantage  of  the  method. 
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Plant 


OUTPUT  PREDICTION  ERROR 


EXTERNAL  INPUT 


INPUT  TO  THE  PLflNT 


APPENDIX  A 


DESCRIPTION  OF  COMPUTER  PROGRAMS 

Computer  programs  are  written  using  WATFIV  and  FORTRAN 
programming  languages.  All  of  them  are  prepared  as  an 
interactive  program.  The  following  explanation  is  about  how 
one  can  use  these  programs. 

The  program  named  CONDIS  given  in  Appendix  B  finds  the 
pulse  transfer  function  from  the  continuous  time  transfer 
function.  The  program  asks  for  the  sampling  interval,  orders 
and  coefficients  of  the  continuous  time  transfer  function 
polynomials.  The  continuous  time  system  is  described  in  state 
variable  form  as 

x  =  Ax  +  Bu 

and  the  corresponding  discrete  time  system  as 

x(k+l)  =  <I>x(k)  +  Tu  (k) 

This  program  gives  the  $  and  r  matrices  corresponding  to  A 
and  B . 

There  is  no  limitation  about  the  order  of  the  system  in 
all  programs.  For  higher  order  systems,  the  dimensions  of 
the  declared  matrices  should  be  increased. 


The  program  named  RCIOP  given  in  Appendix  D  simulates 
an  indirect  adaptive  control  system  using  the  projection 


algorithm.  The  program  asks  some  questions  to  the  user  in 
the  following  order: 

1.  The  order  of  a  plant  (N) ; 

2.  The  constants  a  and  c  that  are  given  in  Equation  (3.38); 

3.  Coefficients  of  the  numerator  polynomial  of  the  plant 
(in  ascending  order  of  z)  ; 

4.  Coefficients  of  the  denominator  polynial  of  the 
plant  (in  ascending  order  of  z) ; 

5.  Coefficients  of  the  controller  polynomials  g(z)  and 
P*(z). 

The  external  input  is  defined  inside  the  program.  It 
can  be  changed  when  it  is  necessary. 

The  last  program  RCIOR  given  in  Appendix  C  is  the  most 
important  one.  Least-squares  algorithms  are  used  for  parameter 
estimation.  This  program  asks  the  same  questions  about  system, 
except  question  #2. 

Computer  programs  RCIOR  and  RCIOP  give  the  graphical 
and  tabulated  results.  Graphical  outputs  are  obtained  using 


DISSPLA. 


APPENDIX  B 


v  v  ?rm . 


COMPUTER  PROGRAM  CONDIS 


WATFIV  ONUK , XRE  F , EXT 
DEFINITION  OF  VARIABLES: 

INTEGER  N.M,R, I,K,L,L1,L2, II, J, 13, Ml, M2, N3 
integer  M3,L3,Io 
REAL  AA< 2,2) ,BBC2 
REAL 

real  _  _  __  _  __  __ 

REAL  C4U0;  10)  I  C5(  10)  10) ) C6(  10)  10)  'Ltl\  10)  10) 


jer  raj  ,  Li  j  ,  j.  o 

AA{  2,2) , BB(  2 ,  1) , ID( 2 ,2 ) , PS{ 10, 10) 

fi  i6,16),gt( i6,io),6i( 1o,io),c2( io,io) 

c3( 10,10), c8( 10,10  ,d4{ 10,10  ,d9{ 10,10 
C4( 10)l0)  C5(  10!  10))C6fl0'  10) )C7(10' 10) 

REAL  C10(l6,l6) ,C11( l0,l6),D2(l0,l6),D3(l0,l0) 
REAL  D5f  10,-10), D6{  10,  l6)  ,62(  lb.  l6)  ,D8(  10.  l6) 
real  cl4(  l6,l6),cl9(  lo, 16) , c20(  10, lo ) , e4[ 10,10) 

- - 10, 10), Dll  10,10  ,C12(  10,10),C13( l6,l6) 

(10. 10).C16(  10, 10 ) ,  C17(  10, 10),C18(10,  10) 

Of  16, 16) .wit  1, io)/io*o./. tDi( io.io) 

10,  lo  j ,  El(  io,  10  );fe2(  10,10  J,E3(  lo): 


REAL  DIO  10,10 
REAL  C15(10 , 10 
REAL  CllO(.l6 
REAL  CC ' ' 


REAL  D1 
REAL  E5 
REAL  FF 
REAL  TC 


- -  ,  , - -  _  , - 0  /  U<b  I  XV/  /  XV  I  .  L1W  I  iW  ,  X 

10. 10), E6( 10,10), C9( 10,10), ss  2, 2 ) 
2,1) , GG( 1 , 1 ) , ALPHAf 10) , CO( lO), aLp , TAR 
TRG,  TR1 ,  TR2  ,  lR3 ,  TR4,  TR5 ,  TR6 ,  T&7 ,  T^8 


ivunu  X  t  X  ,  liVX  ,  X  XVX*  /  X  XV  ^  /  X  AT  /  Xl\w/  /XX 

real  TR9 , TRlO , TR11 . TK1 , c21(_l0,  l6) 

PRINT,  THIS  PROGRAM  FINDS  THE  Z-DOMAIN  PULSE  TRANSFER' 
PRINT  'FUNCTION  WHEN  GIVEN  THE  S-DOMAIN  TRANSFER  FUNC. 
PRINT, '  T 

PRINT, 'THE  S-DOMAIN. TRANSFER  FUNCTION  SHOULD  BE 
IN  THIS  FORM' 

H( S)=NUM( S)/DEN( S) 

NUM(S)=T(M)S**N-1+ . +T(  2  )  S+T(  1 ) 

DEN(  S)=S**N+D(N)S**N-1+ _ +D(2)S+D( 1) 

PRINT,. 

PRINT. 'ENTER  THE  NUMERATOR  POLYNOMIAL  DEGREE' 

READ, Ml  . 

PRINT, ;  ' 

PRINT, 'ENTER  THE  DENOMINATOR  POLYNOMIAL  DEGREE' 

READ, N 

M2=Ml+l 

DO  2  11=1. M2. 

PRINT, ' Tf ' . II, ' )=?' 

READ , TT( I 1 ) 

CONTINUE 

N3=N-1 

M3=Ml+2 

IF  (Ml.  LT.N3)  THEN 


DO  82  I 1-M3 , N 
TT(  1 1  )=0.  0 
CONTINUE 
END  IF 

D°  PRINT, ,11, ' )=?' 

READ , DD( I 1 ) 

CONTINUE 
DO  22  Il=l,N3 
J=1 

AA( II, J)=0.  0 
CONTINUE 
DO  32  J-2 , N 

DO  102  11=1, N3 

IF~( J^EQ.  13 )  THEN 
AA(  11, J)=l. 0 
ELSE  • 

AA(  II, J)=0. 0 
END  IF 
CONTINUE 
CONTINUE 
DO  72  J=1,N 

AA(N, J)=-DD(  j) 
CONTINUE 
DO  52  I1=1,N3 
J=1 

BBtll, J)=0. 0 
CONTINUE 
BB(N,  J)=1,0 
DO  6^  I 1=1 ,N 

CC(  Il)='tT(  II) 
CONTINUE 
R=1 

CALL  FAMILYfN,^ ,'AA' ) 
CALL  RESULTi  N,M, I.K^AA) 
CALL  FAMILY  N,R*BB'} 
CALL  RESULTi;  N,  R,  J1 ,  W , E 
CALL  FAMILYi  R,N, ;CCT) 
CALL  RESULT* R,N, I, K,CC) 
DO  14  J=1,N 
DO  15  t=l,N 

IF  (I.feQ.  J)  THEN 
ID(  I, J)=l.  0 
ELSE  %  „  „ 

ID( I, J)=0. 0 
END  IF 
CONTINUE 

r'riM'TT'MTTB! 


CALL  FAMILY(  N,M,  ID  ) 

CALL  RESULT(  N,M, I , J, ID) 

M=1 

DO  285  1=1  ,N 
DO  286  J=1,N 

RR( I, J)  =  ID(  I, J) 

CONTINUE 
CONTINUE 
TR1=1.  0 

CALL  ALI(  N.N,  I ,  K,  AA,  Wl) 

w  Wii  i  v,i\  Wii  w  Vih. 

PRINT^TC 
THl=TRl/2 
PRINT, TH1 
TH2=TRl/4 
PRINT,  TH2* 

TH3=TRl/8 
PRINT, TH3 
TH4=TR1/16 
PRINT, TH4 

IF  (TR1.  LE.  TC)  THEN 

TR2=TRl/2 

TR3=TRl/3 

TR4=TRl/4 

TR5=TRl/5 

TR6=TRl/6 

TR7=TRl/7 

TR8=TRl/8 

TR9=TRl/9 

TR10=TR1/10 

TR1 1=TR1/1 1 

CALL  SCALARf  N,N, I , K, TR1 , ID, Cl ) 

CALL  SCALARi  N, N, I , K, TR2 , AA,D2 ) 

CALL  CALCULi  N,N,N, I,K,K,C2,C1.D2) 

CALL  SCALARi  N, N, I , K, T&3 , AA,D3 ) 

CALL  CALCULf  N , N , N , I , K , K , C3 , C2 . D3 ) 

CALL  SCALARi  N, N, I , K, TR4, AA, D4) 

CALL  CALCULI  N,N,N, I, K,K,C4,C3.D4) 

CALL  SCALARi  N, N, I , K, TR5 , AA, D5 ) 

CALL  CALCULI  N,N,N, I, K, K, C5 , C4. D5 ) 

CALL  SCALARi  N,N, I , K,TR6, AA,D6) 

CALL  CALCULi  N, N,N, I , K, K, C6, C5 . D6 ) 

CALL  SCALARi  N , N , I , K , TR7 , AA , D7 ) 

CALL  CALCULi  N, N, N, I ,  K,K,  C7 , C6 . D7 ) 

CALL  SCALARi  N, N, I , K, TR8 , AA, D8 ) 

CALL  CALCULi  N,N,N,  I,  K,K,C8,C7,D8) 

CALL  SCALAR  N,N, I, K,TR9,AA,D9) 

CALL  CALCUL(  N,N,N, I,K,K,C9,C8, D9 ) 


10)) 


.v.v3 


■>v 

s.sS 
\’.V> 
v -V 


•.*  v ^  -'.v.v'.-vv.v.v.v.v • 


CALL  SCALAR!  N,N, I , K, TR10.AA, DIO ) 

CALL  CALCULI  N ,  N ,  N , I , K,K, CIO, C9 , DIO ) 
CALL  SCALAR!  N,  N, I , K, TRll , AA, Dll ) 

CALL  CALCUL( N,  N, N,  I , K,K,C11 , CIO, Dll ) 
CALL  SUM(N,N, I,K,C1,C2,C12) 

CALL  SUMi  N,N, I , K, Clz , C3 , C13 
CALL  SUM  N,N, I,K,C13,C4,C14 
CALL  SUM!  N,N, I,K,C14,C5,C15 
CALL  SUM  N,N,I,K,C15,C6,C16 
CALL  SUMi  N, N, I , K, C16, C7 , C17 
CALL  SUMi  N,N, I , K, C17 , C8, C18 
CALL  SUMi  N,N,I,K,Cia,C9,C19 
CALL  SUM  N,N, I , K, Cl 9.C10, CllO) 

CALL  SUMiN,  N,  I ,  K,  CllO,  Cll ,  PS  ) 

CALL  CALcUfj(  ft, ft,  ft,  I ,  K, K, C21 , AA, PS) 
CALL  SUMIN.N, I,K, ID,C2l,FI) 

CALL  FAM I L^ I  ft ,  ft ,  *  PS ’ T 
CALL  RESULTi  N,N,I,K(PS) 

CALL  FAMILY!  N, N,  ;FI  ) 

CALL  RESULT  N,N,I,K,FI) 

CALL  CALCULi  N, R, N, I , K, K, GT, PS, BB) 
CALL  FAMILY  N,R,;GTM 
CALL  RESULT! N, R, I , L,GT) 


ELSE 
TRG=TR1 
L1=0. 0 
tv?—'?  n 

TRG=TR1/(2**L1) 

IF  (TC.  LE.  TRG)  THEN 
L1=L1+1 
GO  TO  677 
END  IF 
PRINT . LI 


T'R=TRi /(  2**Ll ) 

CALL  s6ALAR( N , N , I , K , TR , ID , Cl ) 

CALL* SCALaI{n , N , I , K , TR2 , AA , C2 ) 
TR3=(  TR**3)/6 

CALL  CALCULI  N,N,N,I,K,K,Dl,AA,AA) 
CALL  SCALAR! N,N, I , K, TR3 , D1 , C3 ) 
TR4=( TR**4)/24 

CALL  CALCULI N,N,N, I , K,K,D2, AA,D1) 
CALL  SCALAR!  N, N, I , K, TR4, D2 , C4 ) 
TR5=(  TR**5 ) /l20 

CALL  CALCULI  N , N , N , I , K , K , D3 , AA . D2 ) 
CALL  SCALAR! N.N, I, K,TR5,D3,C5) 
TR6=(  TR**6)/720 

CALL  CALCULI  N , N , N , I , K , K , D4 , AA . D3 ) 
CALL  SCALAR! N. N, I, K,TR6,D4,C6) 
TR7=(  TR**7 )/5u40 


CALL  CALCUL(  N,N,N, I , K,K,D5 , AA.D4) 
CALL  SCALAR  N.N. I, K,TR7,D5,C7) 
TR8=( TR**8)/40320 

CALL  CALCUL? N , N , N , I , K , K, D6 , AA  D5 ) 
CALL  SCALAR?N,N, I ,K, TR8 , D6 , C8 ) 
TR9=(TR**9}/362880  „  M  _ v 

CALL  CALCUL? N,  N, N, I , K,K, D7 , AA.D6) 


CALL  CALCUL(  N,  N, N, I , K,K, D7, AA.D6) 
CALL  SCALAR? N,N, I, K,TR9,D7,C9) 
TR10=?TR**lO)/3 628800 
CALL  CALCUL?  N,  N,  N , I , K,K,D8,AA, D7 ) 
CALL  SCALAR?  N , N ,1 , K , TkIO , D8 , ClO ) 
TRll=(TR**lljy 39916800 
CALL  CALCUL?  N,  N ,  N , I , K, K.D9.AA.D8) 
CALL  SCALAR? N, N. I, K,  TR11.D9 ,  CU ) 
CALL  SUM( N,  N, t,K/Cl,C2,Cl2) 

CALL  SUM?  N,N, I ,K,C12,C3 ,C13 i 
CALL  SUM?  N,  N, I , K, C13, C4, C14 ) 

CALL  SUMi  N,N, I , K, C14, C5 , C15 
CALL  SUM?  N,  N, I,K,C15,C6,C16) 

CALL  SUM?  N,N, I , K, C16, C7, C17 ) 

CALL  SUMi'  N,N,I,K,C17,C8,C18) 

CALL  SUM?N,N, I,K,C18,C9 .C19} 

CALL  SUM?N,N, I , K, C19 , Cl6, C20 ) 

CALL  SUMi  N.N, I, K, C20 . Cll , PS ) 

CALL  FAMILY(N,N, 'Pl'!>  „ 

CALL  RESULT?  N,N, I , K, PS) 

DO  255  16=1, Ll 

CALL  SCALAR? N,N,  I.K.TK2,  ID,  El) 

CALL  FAMILY?  N,N,TEl’i 

CALL  RESULT? N,N,  I,K,E1)  ^ 


CALL  CALCUL  N, N.N, I, K, 
CALL  SUM?N,N, I.K.El.E^ 
CALL  FAMILY?*!,*!, ’ESM, 
CALL  RESULT?  N,N,  I,  RE¬ 
CALL  CALCUL  N,N,N, I, K, 


E3 ) 

CALL  CALCUL?  N ',  N ',  N '  I  \  K ,  K ,  E5 ,  E3  , 
CALL  FAMILY?  N,N/  E5f  ) 

CALL  RESULT? N,N, I, K.E5) 

CALL  COPY?N,N, I,K,P$.E5) 

CALL  FAMILY? N, N, 1  PS'jl 
CALL  RESULT? N, N, I, K, PS) 

TINUE 

CALL  FAMILY? N,N, 'PS' ) 

CALL  RESULT?  N,N, I ,K, PS) 

CALL  CALCULi  N,R,N, I. K,K,GT, PS, 
CALL  FAMILY?  N,  R,  GT' ) 

CALL  RESULT? N,R, I, K,GT)  „ 
CALL  CALCUL? N, N.N, I, K,K,E6,AA, 


CALL  FAMILY? N,N, '£5' ) 
CALL  RESULT? N,N, I,K,E5) 
CALL  COPY? N.N, I,K,P$.E5 
CALL  FAMILY? N,N,  PSTj 


E2 , AA, PS ) 


E5,E3,PS) 


CALL 

CALL 

CONTINUE 

CALL 

CALL 

CALL 


CALL 


wnuvuui  L'i  .  L'i  f  i  / 

SUM? N, N, I , K. ID.E6,FI) 
FAMILY? N, N, *  FI' T 


E6 , AA, PS ) 


CALL  FAMILY? N,N, 'FI' ) 
CALL  RESULT?  N,N,I,K,cI) 


END  IF 

CALL  CALCUL(  N ,  M ,  N , I ,  K ,  K , FF , RR , GT ) 

CALL  CALCULI  M ,  M ,  N , I ,  K ,  K , GG , CC ,  FF ) 

PRINT 

PRINT! ’NUMERATOR  COEFF. S  OF  THE  PULSE  TRANSFER 
FUNCTION1 
PRINT,  T 

PRINT , ’ NUMCOEF  ,CO( N) 

DO  14  *1  L=1 ,  N 
EXECUTE  AR 

ALP=-TAR/L  „  . 

CALL  CALCUL(N,N,N, I, J, J,PP,FI.RR) 


CALL  CALCUL(N,M,N, I, J, J,FF,RR,GT) 

CALL  CALCUL( M,M,N, I, J, J, GG, CC, FF) 

CO( L)=GG(  M,  M) 

ALPHA( L ) =ALP 
1  CONTINUE 

DO  150  L=1,N3 

PRINT, TNUMCOEF' ,CO(L) 
l  CONTINUE, 

print,; 

PRINT  ' 

PRINT! 'DENOMINATOR  COEFF. S  OF  THE  PULSE  TRANSFER 
FUNCTION' 

PRINT. '  r 
DO  151  L=1,N 

PRINT, ' 6ENUMCO  ,ALPHA(  L) 

CONTINUE 

STOP 

REMOTE  BLOCK  AR 
TAR=0  0 

CALL  CALCUL(N,N,N, I , K, K, PP, FI , RR) 

DO  148  1=1, N 
DO  149  $=1,U 

IF(  I.  EO.  J)  THEN 
TAR=TAR+PP(  I, J) 

END  IF 

)  CONTINUE 

i  CONTINUE 

END  BLOCK 
END 

***THI S  SUBROUTINE  READS  THE  MATRICES  FROM  THE  DATA 
FILE. *** 

SUBROUTINE  ENTER( J, D, E, F, G) 

INTEGER  J  ,  D  ,  E  ,  F 
REAL  G(J,D) 


DO  60  E=l, J 

READ  50 . (G(E, F) , F=1 , D) 

50  FORMAT (  10F8.  4) 

60  CONTINUE 
RETURN 
END 

C  ***THIS  SUBROUTINE  PRINTS  THE  ARRAYS  AS  AN  MATRIX 
FORM. *** 

SUBROUTINE  COPY( H, O, P , S , G1 , K1 ) 

INTEGER  H,  0 ,  P ,  S 
REAL  Gl(H.O)  ,K1(H,0) 

DO  131  P=i,H 
DO  132  3=1,0 

Gl(  P,S)=K1(  P,  S) 

132  CONTINUE 

131  CONTINUE 

RETURN 
END 

SUBROUTINE  RESULT(  H, 0, P, S, T) 

INTEGER  H  ,  0  ,  P  ,  S 
REAL  T  (  H  .  0  ) 

DO  80  P  =  1  ,  H 

PRINT  70 ,( T(P , S) , S=1 . 0 ) 

70  F0RMAT(  t6\T2<  10X,  10(3X^12.  6)) 

80  CONTINUE 

RETURN 
END 

C  ***THIS  SUBROUTINE  CALCULATES  THE  MULTIPLICATION  OF 
TWO  MATRICES.*** 

SUBROUTINE  CALCUL(  U. V, Y, Z, X, ZX, 2T, W, Q) 

INTEGER  U^.Y.Z.X'ZX 
REAL  ZT(U,vj/^(6/’5r),Q(Y,V) 

DO  93  Z=1,U 
DO  92  X=1,V 
ZT(  Z.X)  =0.0 
DO  91  ZX=1,Y 

91  coZTiZ6X)^T(Z,X).W(Z,2X).Q(ZX,X) 

92  CONTINUE 

93  CONTINUE 
RETURN 
END 

C  ***THIS  SUBROUTINE  WRITES  THE  MATRIX  NAME  AND  ROW, 
COLUMN  NUMBERS.  ***/ 

SUBROUTINE  F AM I LY(  ROW, COLUMN, NAME ) 

INTEGER  ROW, COLUMN 
CHARACTER* 2  NAME, 

PRINT  94. 'MATRIX' .NAME 

94  FORMAT(Tl.  ' 0T , T2 '  10X, T12 , A6 , T19 , IX, T20, A2 ) 

PRINT  95 , *  ROW  NUMBER.  f . r6w 


95 

96 


42 

41 


44 

43 


46 

45 


FORMATS T1 . ’O'  T2 , 10X.T12 . All . T23 , 8X, T31 , 12) 

PRINT  96 . ’ COLUMN  NUMBER.  . COLUMN 

FORMAT(  Tl , TOr , T2 , 10X, T12 , A14, T26 , 5X, T3 1 , 12) 

RETURN 

END 

***THIS  SUBROUTINE  CALCULATES  THE  SUMMATION  OF  THE 
absolute  values  OF  THE  SAME  COLUMN  ELEMENTS.  *** 
SUBROUTINE  ALI(  1 1 , 12 , HI , H2 , G1 , W2 ) 

INTEGER  II, 12, HI. £2 
REAL  W2(  1, 12) ,Gl( II. 12) 

DO  41  62=1, l2 

DO  42  Hl=l, II 

W2(  1 ,  H2  )  =<V2(  1 ,  H2  )  +ABS(  Gl(  HI ,  H2  ) ) 

CONTINUE 

CONTINUE 

RETURN 

END 

***THIS  SUBROUTINE  MULTIPLIES  THE  MATRIX  BY  SCALAR 
NUMBER.  *** 

SUBROUTINE  SCALAR ( K1/K2,G11,G2,P1,G3,G4) 
INTEGER  K1,K2.G11.G2 
REAL  G3(K1,K2)  .  G4(  K1 ,  K2  ) ,  PI 
DO  43  G2=l,Ki> 

DO  44  G11=1,K1 
G4(G11,G2)=G3(G11,G2)*P1 
CONTINUE 
CONTINUE 
RETURN 
END 

***THIS  SUBROUTINE  CALCULATES  THE  SUM  OF  THE  TWO 
MATRICES. *** 

SUBROUTINE  SUM(K3.K4, 15 , K5 , XI , X2 , X3 ) 

INTEGER  K3,K4.I5.k5 

REAL  X1^K3  ^K4^x|  (  K3  ,  K4 ) ,  X3  (  K3  ,  K4 ) 

DO  46' I 5=1. K3 

X3(  I5,K5)=fcl(  15 , K5 ) +X2(  I5,K5) 

CONTINUE 

CONTINUE 

RETURN 

END 


SENTRY 


APPENDIX  C 


COMPUTER  PROGRAM  RCIOF 


****************************************** 

*  THESIS  PROGRAM  * 

*  INDIRECT  ADAPTIVE  CONTROL  * 

*  PARAMETER  ESTIMATION  * 

*  USING  REC.  LEAST  SQUARE  ALGORITHM  * 

a***************************************** 

’k'k'k'k'k'k'k-k-k’k'k’kic’k'k-k’k'k'k’k'k'kit’k’fe'k'kie'k'kit’kie'kifk'k'k'k'k’kif 

*  DEFINITION  OF  VARIABLES:  * 

•k'k’k-k'kJcick'k'kit’k'k'k'k'k'k'k'klc'k'k'k’k'k'k'k'k'k'k’k'k'k'kic'k'k'k'k'k'k'k 

INTEGER  N1 ,  N, L, K, LI ,  I '  M, L2 , J, N3 , N2 , J1 , J2 , J3 , Ml 
INTEGER  N4, N$ , N6 . COLuft, ^5 , M5 , J7 , J8 , J9 , L5 , M6 , L7 
INTEGER  DEG1 , DEG2 , DEG3 , K9 , M7 , L8 . M2 , M3 . L4 , M9 . M8 


mm 


1111 


REAL  TSM.TSP . TSF, TSN, PI ,V(200), YM(  200) .MU 
REAL  A( 26, 20  j ,H( 20,20) ,  B(  20 ) , C( 20 . 20 ) . 2( 
REAL  D{10)  .T(20)  .F(20i,G(20j,sl(l6J,S2(  10] 
REAL  U3(  4,1)'C5(  10),C6(lO),MK(  50) ,U4( 4, 1), 
REAL  SUM2  .TElfa?.  SIGMA.  ASIGMA, DY, UP(  200),  Oil 
REAL  01(  56)  ,03(  50 )  .  c6(  100)  .  R(  20)  .  KD(  l66]M 
REAL  Sl(  20) , Hl( lOO) .H2( 100) ,H3(  160) ,  H4(  16( 
REAL  Ul(  4, 1 ) , P3(  4, 1  j ,  AB(  4,4), C4(  10),h5(  IOC 
PI=3. 14l5$26$4 
DO  7777  1=1,2 
Hl(  I  )=0.  6 
H2i  I  =0.  0 
H3i  I  i=0.  0 
H4(  I  =0.  0 
H5i  I  =0.  0 
H  CONTINUE 

DO  20  I  =1,200 
DO  8  J  =1,4 
SE( I, J)=0. 0 
S  CONTINUE 

)  CONTINUE 
M=1 

*  THIS  PART  OF  THE  PROGRAM  ASKS  THE  SYSTEM  * 

*  ORDER  ,  PLANT  AND  MODEL  NUMERATOR, DENOM I. * 


>2(4,4) 


C3(  10 ) 
100) 


\  v  /-  -■ 


f  *  s  '  J  * .■  -*  ^  v*  *  ■  V  1  jl' 


2002 


*  POLYNOMIAL  COEFFICIENTS. 
******************************************** 

WRITEf 6.701) 

FORMATf 1 1Y. 'ENTER  THE  ORDER  OF  THE  SYSTEM’ ) 

READ(  5 ,  * )  N 

N1=N+1 

N2=2*N 

N3=N+2 

N6=2*N+ 1 

M7=10*N*2 

L8=10*N2 

•  CALL  FRTCMSf  ’ CLRSCRN  ’ ) 

WRITE? 6.2002) 

2  FORMAT ( 1  Y, T ENTER  PLANT  NUMERATOR  POLYNOMIAL 
COEFFICIENTS7) 


2003 


2004 

2001 


2005 


2006 

2007 


2009 

2008 


WRITEf  6  2003 ) 

FORMATf 1  ’  .  ’  IN  ASCENDING  ORDER  OF  Z’ ) 

DO  2001  1=1, N 
J7=I-1 

WRITEf 6 . 2004)  J7  .  ,  , 

•  FORMAT( 1  Y,TCOEFF.  OF  Z**( ’ , 12 , ’ )=’ ) 
READ(5>)  65(1) 

CONTINUE 

CALL  FRTCMSf  ’ CLRSCRN  ’  ) 

WRITEf 6.2005) 

FORMAT(  ’  Y , T ENTER  PLANT  DENOMINATOR  POLYNOMIAL 
COEFFICIENTS’ ) 

WRITEf 6.2006) 

FORMATf ’  t,tIN  ASCENDING  ORDER  OF  Z’ ) 

WRITEf 6,2007) 

FORMATf 1  r.' HIGHEST  COEFF.  SHOULD  BE  1.  ') 

DO  2008  1=1, N 
J8=I-1 

WRITEf 6,200?)  J8  ... 

FORMATf '  T,r60EFF.  OF  Z**( ’ , 12, ' )=’ ) 

READf  5 ,  *  )  66(1) 

CONTINUE 
DO  112  L=1 , 200 

VP( L)=SIN( L*PI/3)  +  SIN(  L*PI/4)  +  SIN( L*PI/5) 

+  +SIN( L*PI/7)+SIN( L*PI/2)+SIN( L*PI/6) 

CONTINUE 
DO  243  L=1 , 200 
VC( L)=100. 0 
CONTINUE 
DO  124  L=1,N 
Y(L)=0.  6 
U(.Lj=VP(  L) 

CONTINUE 
Y(N1  )=0.  0 
DO  2010  1=1, N 


«  ^iXC8,U*Y( N1 >  "C6(  1 )  *Y< 1 ) +c5< 1 )  *u( 1 ) 

2010  CONTINUE 

DO  2011  1=1, N 
J9=N1-I 
L5=N6- I 

U3(  J9,M)=C6( I) 

U3(  L5,M  =C5( I) 

2011  CONTINUE 

DO  113  I=1,N 
L1=N+ I 
L2=N1-I 
Fl( I,M)=0. 0 
Fl(Ll,M)=VP(L2) 

113  CONTINUE 

DO  121  1=1, N2 
DO  122  K=1,N2 

IF( I. EQ. K)  THEN 
P2(  I,K  =30.  0 
ELSE 

P2(  I ,K)=0. 0 
END  IF 

122  CONTINUE 
121  CONTINUE 

C  CALL  FAMILY(N2,N2,  P2  ) 

C  CALL  RESULT  N2,N2, I. K,P2) 

CALL  COPYM(N2,N2, I,L,AB,P2) 

DO  123  1=1, N2 
Ul( I,M)=0.  0 

123  CONTINUE 

CALL  TRANS(N2,M, I, K, FI, FT) 

CALL  FRTCMS( Y6l&s6rN  T) 

WRITE(6 . 702 ) 

702  FORMATS ’  y,ENTER  Q(Z)  IN  DESCENDING  ORDER  OF  Z' ) 
WRITE(6.703> 

703  FORMAT/  Y /HIGHEST  COEFF.  MUST  BE  1.') 

WRITER  6. 704) 

704  FORMAT/  y, f*REMARK:  Q( Z)  SHOULD  BE  STABLE  POLYNOMIAL') 

DO  157  1=1, N 

WRITE(  6. 705 }N-I  ,  .  , 

705  FORMAT/  Y/COEFF.  OF  Z**(  ’  ,  12,  ’  )  =  ’  ) 

READ{5 ,  *  )  63(1) 

157  CONTINUE 

CALL  FRTCMS( ’ CLRSCRN  1 ) 


WRITE 

FORMA1 


CALL  FRTCMS( ’ CLRSCRN  ’ ) 

WRITE(6 . 706 ) 

FORMAT( *  Y  ENTER  PSTAR(Z)  IN  ASCENDING 
OF  Z  ) 

WRITE(6, 707 ) 

FORMAT/  t  ,  1  HIGHEST  COEFF.  MUST  BE  1.  '  ) 


'  T,  ENTER  PSTAR(Z)  IN  ASCENDING  ORDER 
OF  Z'  ) 


-  * J  *v 

I 


* REMARK- 1: PSTAR( Z)  IS  THE  DENOMINATOR 


- 

».  *  “ 

'.V/. 

<w.y\ 

I-,/ 

•  «■  y 

•i'.s-;/ 

.Oo 


POLYNOMIAL' ) 

WRITEf 6.709) 

FORMATS  T,  'OF  MODEL') 

WRITE(  6. 710) 

FORMAT(  r '*REMARK-2:  PSTAR( Z)  SHOULD  BE  STABLE 
POLYNOMIAL ' ) 

DO  158  1=1, N 

WRITEf 6. 711 )N- I 

FORMAT* *  ^  ,/COEFF.  OF  Z**(  ' , 12 , '  )  = '  ) 

READ* 5  ,  *  )  C4(  I ) 

CONTINUE 


DO  241  1=1.20 
V*  IJ=VP(  I) 
CONTINUE 
DO  114  J=N1,L8 


CALL  CAL6UL( N2 . M , N2 , I , L . K, VI ,  P2  ,  FI ) 
CALL  CALCUL(M.M,N2, I , L, K, V2 , FT,  VI ) 


TSP=1. +V2 
CALL  CAL C 


nijw  wui  i  t  < 

+  V2(-M.M) 
ALCUL(M, 
J)-V4  M, 


M,  N2 , I , L, K, V4, FT, U1 ) 


TSF=Y[j)-V4  M,M) 

TSN=TSF/TSP 

CALL  SCALAR( N2 , M, I , K, TSN, VI , V5 ) 

CALL  SUM*N2,M. I ,K,U1, V5,U21 
CALL  CALCUL* N2 , M, N2 , I , L . K, P3 ,P2, FI ) 
CALL  CALCUL/M ,  M,  N2  ,  I ,  L,  K,  P4,  FT,  P3  ) 
TSM=l./(  1  +  P4f  M,M)7 

CALL  CALCUL* M.N2 , N2 , I , L, K, P5 , FT, P2 ) 
CALL  CALCUL  N2 ,  N2 ,  M,  I ,  L, K . P6, F1.P5} 
CALL  CALCUL* N2 , N2 , N2 . I ,  L, K, P7, P2 , P6 ] 
CALL  SCALAR  N2,N2, I, L,TSM,P7,P8) 
CALL  SUB(N2,N2, I,L,P2,P8,P9) 
INITIALIZE 
DO  115  1=1. N2 
DO  116  K=1,N2 
P2(  I , K)=0. 0 
CONTINUE 
CONTINUE 

CALL  COPYM( N2 , N2 , I , L, P2 , P9 ) 

IF*  J.  EQ.  10)  T6EN 

CALL  C0PYM(N2,N2, I,L,P2,AB) 

ELSE  IF(J.  EQ.  20}  THEN 
CALL  COP YM( N2 ,N2 , I , L , P2 , AB ) 

ELSE  IF*  J. EQ. 30}  THEN 
CALL  COPYM* N2 , N2 , I , L , P2 , AB ) 

ELSE  IF* J. EQ. 40}  TftEN 
CALL  COPYM* N2.N2 , I , L, P2 , AB) 

ELSE  IF*  J. EQ. 50}  THEN 
CALL  COPYM(N2,N2, I, L,F2, AB) 

ELSE  IF* J. EQ. 60}  THEN 
CALL  COPYM* N2 , N2 , I , L , P2 , AB ) 


THEN 


12 ,  N2 ,  I , L, P2 , AB) 
GO  TO  9591 


7011 


4001 


7007 


1T(  N2 ,  N2 , 1 ,  K,  P2  ) 


ELSE  IFfJ.EQ.  70}  THEN 
CALL  COPYM( N2  ,  N2  ,  I ,  L ,  P2  , AB ) 

END  .IF 

I F(  J.  LE.  M7)  GO  TO  9591 
C5( 1 ) =3 . 

DO  2031  1  =  1 ,  N 
J9=N1- I 
L5=N6-I 

U3(  J9 , M )=C6( I ) 

U3(  L5,M  =C5( I 
2031  CONTINUE 
9591  DO  117  1=1 , N2 

Ul( I , M)=0. 0 
117  CONTINUE 

CALL  COPYMf N2 ,M , I , L  ,U1 . U2 ) 

C  CALL  FAMILY(N2,N2, rP2T) 

C  CALL  RESULT  N2,N2, I, K,P2) 

DO  148  1=1,  N 

u(  i )  =v(  i ) 

148  CONTINUE 

DO  211  1=1, N 
D(  I )=Ul( I , M) 

211  CONTINUE 

DO  212  1=1 , N 
L4=N+ I 

B(  IJ=U1(  L4,M) 

212  CONTINUE 

DO  7011  1=1 , N 
M9=N- 1+ 1 
01(  I )=C3( M9 ) 

7011  CONTINUE 

Ql(  N1 )  =1.  0 
$2(  N1 1=0. 0 
DO  4001  1=1, N 

03< I)=D(  t )-C4( I) 

4001  CONTINUE 

DO  7007  1=1 ,N 
M8=N-I+1 
02( I)=Q3(M8) 

7007  CONTINUE 

DEG1=N+ 1 
DEG2=2*N 
DEG3=DEG2+2 
DO  6001  1=1 , DEG3 
R( I}=0. 0 

DO  6002  K=1 , DEG2 

IF(  (_  (  I-K^.  LT.  0.  0).  OR.  (  (  I-K) .  GT.  DEG1 )  ) 
GO  TO  6002 

IF((  I.  GT.  1).  OR.  (  K.  GT.  1}}  GO  TO  6503 
6503  R(  1-1  )=R(  I-1)+Q1(K-1)*02(  I-(K-l)  ) 


AD-A172  351 


ADAPTIVE  CONTROL  MITH  FINITE  TINE  PERSISTENCY  OF 
EXCITATION(U)  NAVAL  POSTGRADUATE  SCHOOL  HONTEREV  CA 
I  ONUK  JUN  86 


UNCLASSIFIED 


F/Q  S/2 


6002 

6001 


9003 


CONTINUE 

CONTINUE 

gSWOT1' 

K9=N6-I+1 
F(  IJ=R(  K9) 

CONTINUE 

N4=2*N+1 

M1=N+1 

DO  41  J4=l , Ml 
DO  51  I=1,N4 
K=I-J4 

A(  I ,  J4)=0.  _ 

IF  (K.  EQ.  0.  )  THEN 
A(  I,  J4)=l. 

END  IF 

IF  ((K.GT.0.  ).  AND.  (  K.  LT.M1))  THEN 
X(I/J4)=D(K) 

END  IF 
CONTINUE 
CONTINUE 
DO  61  J4=1,N 
DO  71  1=1, N4 

A( I , N+1+ J4)=0. 

K=I-J4-1 

IF  ((K.  GT.  0.  ).AND.  (K.  LT.  Ml))  THEN 
A(I  ,N+1+ J4)=B{  K) 

END  IF 
CONTINUE 
CONTINUE 

DO  10  COLUM=l , N4 
SUM2=0. 

DO  70  J4=COLUM,N4 

SUM2=SUM2  +  (  A(  J4, COLUM) ) **2 
CONTINUE 
SIGMA=SORT(  SUM2) 

ASIGMA=XbS(  SIGMA) 

IF  TaSIGMA. LT.  0.  1)  GO  TO  6666 

FORMAT/  ^1^’!  'MATRIX  A  IS  SINGULAR') 
END  IF 

DO  80  J4=C0LUM,N4 
G(  J4J=0. 

CONTINUE 

G( COLUM ) =S I GMA 

SuM2=0. 

DO  90  J4=COLUM,N4 

SUM2=SUM2+( A( J4, COLUM )-G( J4) )**2 
CONTINUE 

DO  100  J4=C0LUM, N4 


6666 


DO  110  K=C0LUM,N4 

IF  ( ABS( SUMi).  LT. 0.  00001)  THEN 
TEMP=0. 

ELSE 

TEMP=-2*(  A(  J4, COLUM) -G( J4) )*( A( K,COLUM) 
►  -G( K) )/s6M2 

END  IF 

IF  (J4.  EQ.  K)  THEN 
TEMP=I.  +TEMP 
END  IF 

H( J4,K)=TEMP 
CONTINUE 
CONTINUE 

DO  131  J4=C0LUM, N4 
DO  130  K=COLUM,N4 
TEMP=0. 

DO  140  L=C0LUM,N4 

TEMP=TEMP+H(  J4,L)*A(  L,K) 

CONTINUE 
C( J4,K)=TEMP 
CONTINUE  • 

CONTINUE 

DO  150  J4=C0LUM,N4 
TEMP=0. 

DO  160  L=COLUM.N4 

TEMP=TEMP+H(  J4,L)*F(  L) 

CONTINUE 
T(  J4)=TEMP 
CONTINUE 

DO  170  J4=C0LUM, N4 
DO  180  K=COLUM , N4 

co^feK,=C(J4'K) 

F(J4)=T( J4) 

CONTINUE 

CONTINUE 

X(  N4)=F(  N4)/A(  N4,N4) 

I— N4-1 
SUM2=0. 

M2=I+1 

DO  30  J4=M2.N4 

SUM2=SUM2+A( I, J4)*X( J4) 

CONTINUE  „  %  . 

1 11 1  A^Ii  il  1  J^L^o!  060&1 )  THEN 
X{  I)=0. 

END  IF 

IF( I. GE. 1)  GO  TO  812 
Hl( J)=X( 1) 


A  ;  V  .•  V.  /.V.  ,'. 


H2(  J )=X( 2 
H3( J  =X  3 


H4{  J)=X{  4) 

H5  J)=xl 5} 

DO  141  1=1. Ml 

con?^=x(i) 

DO  142  I =N3 , N4 
N5=I-N1 
Sl(  N5)=X(  I) 

contiMje 

UP( J)=0. 0 
DO  5001  1=1 ,  N 

UP(J)=UP^  J W  S2(L7  j  -C3(  I )  )  *U(  J-I  )+Sl(  I )  *Y(  J-I ) 
CONTINUE  + 

US5=jllP( J>+V( J) )/( 1"S2( 1} ) 

Y( J5 )=0. 0 
DO  3603  1=1  ,N 
M6=N1- I 

Y(J5)=Y(  J5  )  -C6(  I )  *Y(  J5-M6  )+C5(  I)*U(  J5-M6) 
CONTINUE 
DO  189  1=1 ,N1 
YM(  I  )=0.  0 
CONTINUE 
YM( J5 )=0.  0 
DO  2508  1=1  ,N 
M8=N-I+1 
Bl(  I1=C5(M8) 

CONTINUE 
DO  3004  1=1. N 

CONTliui)™^  J5)-C4(  I)*YM(  J5-I  )+Bl(  I  )*V(  J5-I) 

MK( J5)=-Y(  J5 J  - 
DO  3005  1=1.  N 

MK(J5  )  =Mfc(  J5  )+B(  I )  *  V(  J5- 1 )  -C4(  I )  *Y(  J5  - 1 ) 
CONTINUE 

MU( J5 ) =ABS{  MK(  J5 ) ) 

IF(MU(  J5).  LT.  20.  0)  THEN 

el^)=^C(J5) 

V(  J5)=VC(  J5 )+VP(  J5 ) 

END  IF 

DO  120  1=1 ,N2 
Fl(  I . M)=0.  0 
CONTINUE 
DO  118  1=1, N 
L1=N+ I 


3001 


3002 


8888 

C 

C 

Clll 

C169 

1111 


8002 

976 

C 

C 

C 

C 

C 

c 

4445 

4442 

4441 


1025 

C  **** 


Fl]  Li y  A)=u{  J2 
CONTINUE 
DO  119  1=1, N2 
FT( M, I)=0. 0 

contiMje 

CALL  TRANS( N2,M, I . K, FI, FT) 
CALL  FAMILY(  N2,M,  ^Ul' ) 

CALL  RESULT(  N2,M,  I,K,U1) 

do  3001  i=i;n2 

U4(  I,M)=0.  0 
CONTINUE 

CALL  SUB(N2,N2,  I ,  L,  U3  ,  Ul, U4) 
C7(J)=0.  0 
DO  3002  1=1. N2 

C0N?M=C7(J)+U4<I'M>*U4(I'M) 

g§^isggRT(c7(J>) 

DO  8888  1=1, N1 
C8( I )=0. 

CONTINUE 
DO  8111  1=1.100 
WRITE  (8,8169)  C8(I) 

CONTINUE 


FORMATS 10X.F12.  6) 

WRITEfb. llllT'SYs.  OUTPUT' , 'MOD.  OUTPUT' 
FORMAT(  i3X, All , 6X, All) 

DO  976  1=1.100 

WRITE  (6,8002)  Y(I),YM(I) 

WRITE  (8.8002)  Y(.I),YM(I) 

FORMAT!  9&,  FIS.o,  oX, £l5.  6) 

CONTINUE 
DO  4441  1=1 ,L8 
WRITE(  8,4445) I 
WRITE!  8, 4442  H1(I) 

WRITE!  8,4442  H2  I  i 
WRITE!  8-,  4442'  H3  I 
WRITE! 8,4442  H4(  I 

WRITE! 8,4442  H5  I 
FORMAT( T  ,12) 

FORMAT!  10X,F15. 6) 

CONTINUE 
DO  1025  1=1. L8 

CONTINUE 

************************************** 


THIS  PART  PLOTS  THE  NECESSARY  DATA  * 
THAT  PROGRAM  CALCULATES  * 

If**************************************** 

CALL  VPLOT(KD,SE.L8,2. 'TIME  '  ' 

CALL  VIPLOTf KD,C8, L8, I, ;TIME  , 

CALL  V2PLOTi' KD,MU,  L8.1.  TIME  ' 
CALL  V3PLOTi  KD, V,L8, 1, ’TIME  ' 

CALL  V4PLOT.  KD,U.L8.1.  'TIME  '  ' 

CALL  V4PLOT  KD ^1^8,1,! 

CALL  V4PLOTi  KD ,  H2 ,  L8 , 1  ' 

CALL  V4PLOTi  KD,H3,L8,1  ' 

CALL  V4PLOTi' KD,H4#  L8, 1,  ' 

CALL  V4PLOT( KD,H5, L8, 1, 

STOP 

END 

***************************************** 
THIS  SUBROUTINE  PRINTS  THE  ARRAYS  AS  AN* 
MATRIX  FORM  -  * 

t**************************************** 

SUBROUTINE  RESULT(  H, 0, P, S,  T) 

INTEGER  H  ,  O  ,  P  ,  S 


REAL  T  (  H  .  0  ) 
DO  80  P  =  1  ,  H 


C  ** 
C  * 

C  * 

C  ** 


CONTINUE 

RETURN 

END 

***************************************** 
THIS  SUB.  CALCULATES  THE  MULTIPLICATION* 
TWO  MATRICES.  * 

lr  **************************  ************** 

SUBROUTINE  CALCUL( U . V , Y, Z , X , ZX , ZT , W , Q ) 
INTEGER  U,V.Y,Z.X'Z& 

REAL  ZT(  U,  V  j ,  &(  6,  Y ) ,  Q(  Y,  V) 

DO  93  Z=1  U 
DO  92  X=1,V 
ZT(Z.X)  =0.0 
DO  91  ZX=1, Y 

cont?n6e)=4t( Z,x)<w(  Z- ZX,*2< ZX,X) 

CONTINUE 

CONTINUE 

RETURN 

END 

k'k'k'k'k'k'kick’k'kit'kic'kle'it'k'kie’kic'k'k'k'k'k'k'k'k'k'k’kic'k'k'k'k'k'k'k 

THIS  SUBROUTINE  WRITES  THE  MATRIX  NAME  * 
AND  ROW. COLUMN  NUMBERS.  * 


SUBROUTINE  FAMILY(  ROW, COLUMN, NAME ) 

INTEGER  ROW, COLUMN 

CHARACTER* 2  NAME 

WRITE(  8, 94ji  1  MATRIX  ,  NAME 

PRINT  94, r MATRIX1 .NAME 

PflDMST/  'P'i  1  ftV  T1  O  SC  Q  IV 


FORMATf Tl, ’0; -T2, 10X.T12.A6.T19, 1X,T20,A2) 
WRITE( 8,9$)  ROW  NUMBER. ’ , ROW 
PRINT  95 , 'ROW  NUMBER.  T, ROW 

FORMAT(Tl , 1 0 ' , T2 , 10X, T12 , All , T23 , 8X, T31 , 12 ) 
WRITE( 8,96)  ' COLUMN  NUMBER.  T. COLUMN 
PRINT  96. tCOLUMN  NUMBER.  T. COLUMN 
FORMAT( Tl, TOr , T2 , 10X, T12 , A14/T26, 5X, T31 , 12) 
RETURN 
END 

******************************************** 

*  THIS  SUBROUTINE  CALCULATES  THE  SUMMATION  * 

*  OF  THE  SAME  COLUMN  ELEMENTS.  * 

*  ******************  ******* ★★**★******★★★**** 

SUBROUTINE  ALIfl 1 , 12 , HI , H2 , G1 , W2 ) 

INTEGER  II, 12, HI, H2 
REAL  W2(  1, t2),Gl(  II. 12) 

DO  41  62=1. l2 

DO  42  6l=l.Il 

_ _ „  W2  (  1 , H2 ) =62  (  1 , H2 ) +ABS( Gl( HI , H2 ) ) 

CONTINUE 

CONTINUE 

RETURN 

END 

****************************************** 

*  THIS  SUBROUTINE  MULTIPLIES  THE  MATRIX  * 

*  BY  SCALAR  NUMBLR.  * 

*★******★*★*★★  v- 

SUBROUTINE  SCALAR( K1 , K2 , G1 1 , G2 , PI , G3 , G4 ) 
INTEGER  K1,K2.G11.G2 
REAL  G3(Kl,K2j  . G4(  K1 , K2  ) , PI 
DO  43  G2=1,k2 

DO  44  Gll*l, K1 
G4( Gil , G2 )=63(  Gil , G2 ) *P1 
CONTINUE 
CONTINUE 
RETURN 
END 

****************************************** 

*  THIS  SUBROUTINE  CALCULATES  THE  SUM  OF  * 

*  THE  TWO  MATRICES.  * 

INE  SUM{  K3 ,  K4,  15,  K5 , XI ,  X2 , X3 
K3 . K4 , 1 5 ,  K5 


46 

45 


C 

C 

C 


46 

45 


C 

C 

C 

C 


46 

45 


C 

C 

G 

C 


46 

45 


C 

C 

C 


DO  46  I5=1.K3 

X3( I5,K5)=X1( 15 , K5 )+X2(  15, K5) 

CONTINUE 

CONTINUE 

RETURN 

END 

****************************************** 

*  THIS  SUBROUTINE  COPPIES  THE  MATRICES.  * 
****************************************** 

SUBROUTINE  COPYM( K3 , K4, 15 , K5 , XI ,X2 ) 

INTEGER  K3 , K4, 15,  K5 
REAL  X1^K3  ^K4^X2  ( K3 , K4 ) 

^  DO  46' 15=1. K3 

,mWTWtr  Xl(  I5,K5)=&2(  15, K5) 

CONTINUE 

CONTINUE 

RETURN 

END 

****************************************** 

*  THIS  SUBROUTINE  FINDS  THE  TRANSPOSE  * 

*  OF  THE  MATRIX.  * 

****************************************** 

SUBROUTINE  TRANS( K3 , K4, 15 , K5 , XI , X2 ) 

INTEGER  K3/K4.I5'K5 
REAL  X1(K3,K4J  ,X2(K4,K3) 

DO  45  K$=1,K4 

DO  46  15=1, K3 
X2(K5, 15)=X1(  15, K5) 

CONTINUE 

CONTINUE 

RETURN 

END 

************* *  * ******************************* 

*  THIS  SUBROUTINE  CALCULATES  THE  SUBTRACTION.  * 

*  OF  THE  MATRICES.  -  * 

********************************************** 

SUBROUTINE  SUB( K3 . K4, 15, K5,X1,X2 ,X3 ) 

INTEGER  K3,K4.I5.fc5 

REAL  X1^K3  ^K4 S  X2  (  K3 , K4 ) , X3  (  K3 , K4 ) 

DO  46'l5=l.K3 

X3( I5,K5)=X1( 15, K5 ) -X2(  15, K5) 

CONTINUE 

CONTINUE 

RETURN 

END 

•k'k'k'k'klc’k'k'k'k'k'k'k'k’k’kit'k'k'k'k’k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

*  THIS  SUB.  PLOTS  THE  ARRAYS  USING  * 

*  THE  DISSPLA  PROGRAM.  * 


V-  *rJ 


■  V  *.■  V 

•  V  .*  v 
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m 

»  * 
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V.'v'-y 
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****************************************** 


SUBROUTINE  VPLOT( T,U,M,N, LABELX, LABELY) 

C************* 

REAL  U(800, 4)  ,T(800 ) , Y(  800 ) , YX( 4), YN( 4) 
INTEGER  PMAX, PMlN,M,N, LABELX, LABELY 

CALL  AMAX( T , M , TMAX , PMAX ) 

CALL  AMINf T,M,TMIN,PMIN) 

- J  =l,ft 


CALL  AMINf  T,M,TMIN,PMIN) 

DO  20  J  =l,fc 
DO  10  I  =1 . M 
Yf  I )  =U(l,J) 

CONTINUE 

CALL  AMAXf  Y , M , YMAX , PMAX ) 

YXf J ) ■ =  YMAX 

CALL  AMIN(Y,M, YMIN, PMIN) 

YNfJ)  =YMIN 
WRlTE( 6,11)  YMAX, YMIN 
CONTINUE 

WRITE(6,13)(YX(  I)  ,1=1,4) 

WRITEf  6, 13  J  i  YN(  IJ,  1=1,4) 

CALL  AmAx(YX,N, YMAX, PMAX) 

CALL  AMINf YN.N, YMIN .PMIN ) 

WRITE(6. 11)  YMAX, YMIN 
CALL  TEK  618 
CALL  BLOWUP( 1.2) 

CALL  PAGEf  11.  ,8.  5) 

CALL  NOBRDR 
CALL  AREA2D(  9.  ,6.  ) 

CALL  XNAME( LABELX, 6) 

CALL  YNAMEf  LABELY, 6 ) 

CALL  HEADIN  (' PLANT  AND  MODEL  OUTPUTS  $’,100,1. ,1) 
CALL  CROSS 

CALL  GRAFf TMIN, 'SCALE' , TMAX, YMIN, 'SCALE' , YMAX) 

CALL  SPLINE 
DO  40  J  =1,N 
DO  30  I  =1.M 
^  Y(  I)  =U(!,J) 

CONTINUE 

IFfJ.  EO.  2)  CALL  CHNDOT 
CALL  CURVE  (T,Y,M,0) 

CONTINUE 
CALL  ENDPL(O) 

CALL  DONEPL 

FORMAT (  ,  ] , 10X,2G12.  5/) 

FORMAT('  ' ,10X,4G12.  5/) 

RETURN 

END 


f -r w 

■V/  * 

a^v*- 


as* 


,  10X, 4G12. 5 


SUBROUTINE  V1PL0T( T, U, M, N, LABELX , LABELY ) 


vvo 


*  »  ■  * 


vwi 


vsy 


%•  < 


VSV'i  .■'•VTl.li  kinWAfl  /V  u- 
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20 


30 


40 


CALL  AMAX( T , M , TMAX , PMAX ) 
CALL  AMIN?  T, M, TMIN, PMIN) 

DO  20  J  =1,N 
DO  10  I  =1 , M 
Y(IJ  =U(I,J) 

CONTINUE 

CALL  AMAX(  Y , M , YMAX , PMAX ) 
YX£J)  =  YMAX 


CALL  AMIN(Y,M, YMIN, PMIN) 

YN(  J)  =YMIN 
WRITE( 6,11)  YMAX, YMIN 
CONTINUE 
WRITE(  6 , 13 )  (  YX( 

WRITE?  6,13)? YN? 

CALL  AMAX(YX,N, 

CALL  AMIN? YN,N, YMIN, PMIN 
WRITE(6. 11)  YMAX, YMIN 
CALL  TEK  618 
CALL  BLOWUP (1.  2 ) 

CALL  PAGE( 11.  ,8.  5) 

CALL  NOBRDR 
CALL  AREA2D( 9 .  , 6.  ) 

CALL  XNAME( LABELX, 6) 

CALL  YNAME?  LABELY ,6 j 
CALL  HEAD IN  ('OUTPUT 
CALL  CROSS  .  . 

CALL  GRAF(  TMIN, ’ SCALE' , TMAX, YMIN, ' SCALE' 
CALL  SPLINE 
DO  40  J  =1 , N 
DO  30  I  =1.M 

S^=U(I-J) 

CALL  CHNDOT 
( T, Y,M, 0) 


PREDICTION  ERROR  $' 


11 

13 


CONT 
IF(  J.  EO.  2) 
CALL  CURVE 
CONTINUE 
CALL  ENDPL(O) 
CALL  DONEPL 

UADMAfP/  »  1 


FORMAT ( . 
FORMAT?  ' 
RETURN 
END 


, 10X, 2G12. 

'  , 10X, 4G12. 5 


SUBROUTINE  V3PLOT(  T, U, M, N, LABELX, LABELY) 

Q* ************ 


REAL  U( 800 . 4 ) , T(  800 ) , Y( 800) . YX( 4 ) . YN( 
INTEGER  PMAX, PMIN, M,N, LABELX, LABELY 

CALL  AMAX(T,M, TMAX, PMAX 
CALL  AMIN? T,M, TMIN, PMIN 


4) 


100,1. ,1) 

YMAX) 


105 


>  / 


oooooooo 


CALL  AMAXJ  Y,  M ,  YMAX ,  PMAX ) 

YX( J )  =  YMAX 

CALL  AMIN(Y,M, YMIN,PMIN) 

YN(J)  =YMIN 

WRITE( 6.11)  YMAX,  YMIN 


WRITE( 6,11)  YMAX, YMIN 
CONTINUE 

WRITE(  6 , 13 )  (  YX(  I), 1=1, 4) 

WRITEf  6, 13j  (  YN  IJ  ,  1=1 , 4j 
CALL  AMAX( YX ,  N , YMAX , PMAX ) 

CALL  AMIN( YN, N, YMIN, PMIN) 

WRITE( 6 . 11 )  YMAX, YMIN 
CALL  TEK  618 
CALL  BLOWUP( 1.2) 

CALL  PAGE( 11.  ,8.  5) 

CALL  NOBRDR 
CALL  AREA2D(  9.  ,6.  ) 

CALL  XNAME( LABELX, 6) 

CALL  YNAME( LABELY, 6 

CALL  HEADIN  ( T INPUT  TO  THE  PLANT  $',100,1. ,1) 
CALL  CROSS  , 

CALL  GRAF( TMIN, ' SCALE' , TMAX, YMIN, 'SCALE' , YMAX) 
CALL  SPLINE 
DO  40  J  =1 , N 
DO  30  I  =1 . M 
Y(  I )  =U(I,J) 

CONTINUE 

IF(J.  EQ.  2)  CALL  CHNDOT 
CALL  CURVE  (T,Y,M,0) 

CONTINUE 
CALL  ENDPL(O) 

CALL  DONEPL 

FORMAT(  '  ; ,10X,2G12.  5/) 

FORMAT( '  ' , 10X, 4G12. 5/) 

RETURN 

END 

SUBROUTINE  AMAX(  Y, N, YMAX, PMAX) 

******************************************** 

*  SUBROUTINE  TO  COMPUTE  THE  LARGE  ELELEMENT 

*  IN  A  GIVEN  ROW  OR  COLUMN  IN  ARRAY  TA. 
★*********★*★*★***★***★★****★***★*********** 

***  VARIABLE  DECLARATION  *** 

REAL  AMAX 
INTEGER  PMAX 
DIMENSION  Y(N) 

YMAX=Y(  1) 

DO  5100  1=2. N 


s 


oooooooo  on 


IF(  Y{ I).  LE.YMAX)  GO  TO  5000 
YMAX=Y( I ) 

PMAX=I 

5000  CONTINUE 

5100  CONTINUE 

C  MAXR0W=P0S 

RETURN 
END 


SUBROUTINE  AMIN(  Y, N, YMIN, PMIN) 

******************************************** 

*  SUBROUTINE  TO  COMPUTE  THE  SMALL  ELEMENT 

*  IN  A  GIVEN  ROW  OR  COLUMN  IN  ARRAY  TAT. 
******************************************** 


C 


6000 

6100 


***  VARIABLE  DECLARATION  *** 
REAL  AMIN 
INTEGER  PMIN 
DIMENSION  Y(N) 


YMIN=Y(1] 

DO  6100  1=2 ,N 


IF( Y(I ) .  GE.  YMIN)  GO  TO  6000 
YMIN=Y( I ) 

PMIN=I 

CONTINUE 

CONTINUE 

RETURN 


APPENDIX  D 


COMPUTER  PROGRAM  RCIOP 


C 

c 

c 

c 

c 

c 

c 


****************************************** 

*  THESIS  PROGRAM  * 

*  INDIRECT  ADAPTIVE  CONTROL  * 

*  PARAMETER  ESTIMATION  * 

*  USING  PROJECTION  ALGORITHM  * 

*************************  ****★★★★★★*★★★**★ 

DEFINITION  OF  VARIABLES: 

INTEGER  N1 , N , L - K, LI , I,  M,  L2 , J , N3 , N2 , J1 , J2 , J3 , Ml , M2 
INTEGER  N4 ,  N$ ,  &6 .  COLUM ,  J5 ,  M5 ,  J7 ,  J8 ,  J9 ,  L5 ,  M6 ,  L6 ,  L7 
INTEGER  DE6l , f>EG2 , DEG3 , K9 , M7 , L8 , M9 . M8 . M3 , L4 

m  ajiiihB&iWiMjjkrt  is  kfi'.ij ^2(4,4, 

REAL  ABU,4i,P3(  4,  4},  P4r4.  ij i,  Eli  20}  ,  Ul(4.  1 ),  R{  26 } 
REAL  TSP .  tSF,  TSN .  £  I '  </(  800  S.  YAfSOO )  .  &U(  8o6}  .  fa>{  800 ) 
REAL  A(26,20)  . H( 2o . 2o) , B( 26) , C( 20. 20} . X(  20} 

REAL  D(10} , T( 2o) .  F(  20 ) , G(20Ji , Sl(  l6 ) , S2( 10 ) ,  C3(  10} 
REAL  U3(  4 '  i ) ,  C5(  to}  ,  Cbi  10} ,  Mfc(  800}  ' <74{  4, 1 ) ,  C7(  800 ) 
REAL  VCj86ol.KD(.800Lc8(806),C4(l6hQ3(6o) 

REAL  SUM2/TEm>/  SIGMA ,ASiGMA,DY,uP(  660}. Q2(  50) 

REAL  CONI .C0N2 , se( 806 , 4 ) , v6(  1 , 1 ) , ql(  50 } 

¥1=3. 141592654 
DO  20  I  =1,200 
DO  8  J  =1,4 
SE(  I , J)=0. 0 
CONTINUE 
CONTINUE 
M=1 

WRITE(6,701) 

FORMAT(  1 1Y, 'ENTER  THE  ORDER  OF  THE  SYSTEM’ ) 

«'*>  * 

N2=2*N 
N3=N+2 
N6=2*N+1 
L8=10*N2 
WRITE( 6.2301) 

FORMAT( ’  VENTER  CONSTANT  C  (PROJECTION  ALGORITM)  ) 
READ(5  *}  CONI 
WRITE(  6.2302) 

FORMAT( *  t,tENTER  CONSTANT  A  (PROJECTION  ALGORITHM)’) 
READ(5  *}  CON2 
WRITE(  6,2002) 

2002  FORMAT* ’l  VENTER  PLANT  NUMERATOR  POLYNOMIAL 


8 

20 


701 


2301 


2302 


2003 


2004 
2001 

2005 

2006 
2007 


2009 

2008 


COEFFICIENTS' ) 

WRITEf 6.2003) 

FORMAT*  1  VIN  ASCENDING  ORDER  OF  Z' ) 

DO  2001  1=1,  N 
J7=I-1 

WRITE* 6.2004)  J7  ,  ,  , 

FORMAT* *  r , 'COEFF.  OF  Z**( ',12,'  =') 

READ* 5 ,  *  )  65*1) 

CONTINUE 
WRITE* 6.2005) 

FORMAT* ' 1*T ENTER  PLANT  DENOMINATOR  POLYNOMIAL 

coefficients”) 


FORMAT* '  ,  H 

DO  2008  1=1 ,N 
J8=I-1 


WRITE* 6.2009)  J8 
FORMAT* 1  t7tCOEFF.  OF  Z* 


*(  M2,  '  )  =  '  ) 


2010 


2011 


J-  -fa  I  g  WVSAJfa*  fa  •  Vfa  tJ  V  /AM/  I  / 

READ*  5  ,*)  66(1) 

CONTINUE 
DO  112  L=l, 200 

VP*  L)=SIN(  £,*PI/3)  +  SIN(  L*PI/4)+SIN*  L*PI/5) 
+  +sin(  l*pi/6)  +  sin(  l*pl/7)  +  SIN(  L*PI/2) 
CONTINUE 
DO  243  L=1 ,200 

vc*  lj=i6o.  0 

CONTINUE 
DO  124  L=1.N 
Y(L)=0.  6 
^  uIlUvp*  L) 

CONTINUE 
Y* N1)=0.  0 
DO  2010  1=1. N 

~  61  >  “C6<  I)*V<  I)  +C5*  I )  *U(  I ) 

CONTINUE 
DO  2011  1=1, N 
J9=N1- I 
L5=N6- I 

U3( J9,M'=C6* I) 

U3L5,M)=C5U) 

CONTINUE 
DO  113  1=1, N 
L1=N+I 
L2=N1- I 
Fl( I,M)=0.  0 
Fl(  L1,M)=VP( L2 ) 

CONTINUE 
DO  121  1=1, N2 


DO  122  K=1,N2 

IF(I.EO.  K)  THEN 
P2(  f, K)  =1. 0 
ELSE 

P2(  I ,K)=0. 0 
END  IF 
CONTINUE 
CONTINUE 
DO  1121  1=1, N2 
DO  1122  K=1.N2 

IF(  I.  EO.  ft)  THEN 
P3(  I , K ) =CON2 
ELSE 


1122 

1121 

C 

C 


P2( I ,K)=0. 0 
END  IF 


CONTINUE 

CONTINUE 


DO  123  1=1 ,N2 
Ul(  I  ,M )  =0 .  0 
CONTINUE 

CALL  TRANS(N2 , M, I , K, FI , FT ) 

WRITER 6, 702 ) 

FORMATf’  T. 'ENTER  Q(Z)  IN  DESCENDING  ORDER  OF  Z' ) 
WRITE(6.70$) 

FORMATf'  T, 'HIGHEST  COEFF.  MUST  BE  1.') 

WRITE( 6 . 704) 

FORMATf '  T,'* REMARK:  Q(  Z)  SHOULD  BE  STABLE  POLYNOMIAL') 
DO  157  t=l , N 

WRITEjfS.  705  )N- 1 

FORMAT ( '  r, 'COEFF.  OF  Z**( ' , 12, ' )=' ) 

■  READ(5 ,  *  )  63(1) 

CONTINUE 
WRITE(6 . 706 ) 

FORMAT(  ‘  T, 'ENTER  PSTAR(Z)  IN  ASCENDING  ORDER 
OF  6'  ) 

WRITEjf  6 .707) 

FORMATf'  r. 'HIGHEST  COEFF.  MUST  BE  1.  ' ) 

WRITEf 6 . 706 ) 

FORMAT (  '  T, '*REMARK-1:  PSTAR(Z)  IS  THE  DENOMINATOR 
POLYNOMIAL1 ) 

WRITEf 6.709) 

FORMATf ’  T , ' OF  MODEL ’ ) 

WRITEf 6. 710) 

FORMAT(  '  T. '*REMARK-2:PSTAR(Z)  SHOULD  BE  STABLE 
POLYNOMIAL1 ) 

DO  158  1=1 , N 

WRITE(  6, 711)N-I 


...  -  * 


rv  s 

W-' 

,  V  V 


» *  •  *  » 
■-•tv; 
*  ,*.  * 


VA 

%v-.\ 


*  »  ji  , 

/•A 

w 

v*V*J 


FORMAT*  '  '  'COEFF.  OF  Z**( ' , 12, ' )=' ) 
READ* 5 ,  *  )  64(1) 

CONTINUE 
DO  241  1=1,20 

contS:  NiiiVP  1  * 

DO  114  J=N1,L8 

CALL  CAL6UL( N2 . M,  N2 , 1 . L . K, VI , P2 . FI ) 

CALL  CALCULi  M,  M,  N2 ,  I ,  L,  K,  V2  ,  FT,  VI ) 
TSP=C0N1+V2i  M,m} 

CALL  CALCUL  M, M , N2 , I , L, K , V4 , FT , U1 ) 
TSF=Y(J)-V4  M  M} 

TSN=TSF/TSP 

CALL  CALCUL( N2,M,N2, I , L, K, P4. FI ,  P3 ) 

CALL  SCALAR  N2,M. I, K,tSN.P4,V5) 

CALL  SUM(N2,M,  I  ,K,  U1,V5,U2) 

DO  117  1=1, N2 
Ul( I ,M)=0. 0 
CONTINUE 

CALL  C0PYM(N2.M, I,L,U1,U2) 

CALL  FAMILY(N2,N2,TP2Tj 
CALL  RESULT* N2 , N2 , I , K, P2 ) 

DO  148  1=1, N 
U(  IJ=V(  t ) 

CONTINUE 
DO  211  1=1, N 
D( I )=Ul( I , M) 

CONTINUE 
DO  212  1=1 , N 
L4=N+ I 

B* I )=U1(  L4,M) 

CONTINUE 
DO  7011  1=1 , N 
M9=N-I+1 


01*  I)=C3(M9) 

7011  CONTINUE 

Q1*N1)=1.  0 
02(Nli=0.  0 
DO  4001  1=1, N 

03 ( I )=D( I)-C4( I) 
4001  CONTINUE 

DO  7007  1=1, N 
M8=N-I+1 
02 ( I)=Q3(M8) 

7007  CONTINUE 

DEG1=N+1 
DEG2=2*N 
DEG3=DEG2+2 
DO  6001  1=1 ,DEG3 
R(  I  )=0.  0 


6503 

6002 

6001 


9003 


DO  6002  K=1,DEG2 

IF(  (  io_T^  6002  ‘  °>  •  0R-  (  (  r-K)  •  GT*  DEG1 )  ) 
IF(  (  I.  GT.  1).  OR.  (  K.  GT.  1))  GO  TO  6503 

CONTiNUE)=R(I‘1)+Qi(K"1)*^(I'(K“1)) 

CONTINUE 

WtS&'llVSk1' 

K9=N6-I+i 

coNTiNiiR(  K9) 

N4=2*N+1 

M1=N+1 

DO  41  J4=1,M1 
DO  51  1=1, N4 
K=I- J4 
A(  I  .  J4->=0. 
li? 

SND  IF 


IF  (  (  K.  GT.  0.  )  .  AND.  (  K.  LT.  Ml ) )  THEN 
A( I , J4)=D( K) 

END  IF  ' 

CONTINUE 
CONTINUE 
DO  61  J4=1,N 
DO  71  1=1. N4 

A(_I ,  N+1+ J4)=0. 

™en 

END  IF 
CONTINUE 
CONTINUE 

DO  10  C°LUM=1/N4 
SUM2=0. 

DO  70  J4=COLUM,N4 

CONTINUE™2^  A(  J4'C0LUM>  >**2 
SIGMA=SORT( SUM2 ) 

ASIGHA=ABS(  SIGMA) 

IF  {  ASIGMA.  LT.  0.  00001)  THEN 
WRITE(  6 . 714) 

END^IT^  'MATRIX  A  IS  SINGULAR') 

DO  80  J4=C0LUM. N4 
G(J4J=0. 

CONTINUE 
G( COLUM ) =S I GMA 


J 


SUM2=0. 

DO  90  J4=COLUM,N4 

SUM2=SUM2+( A(  J4, COLUM) -G(  J4) )**2 
CONTINUE 

DO  100  J4=C0LUM,N4 
DO  110  K=C0LUM,N4 

IF  (ABS(  SUM2).  LT.  0.  00001)  THEN 
TEMP=0. 

ELSE 

TEMP=-2*(  A(  J4 , COLUM) -G(  J4) )*( A(  K, COLUM) 
-G(  K)  )/s{jM2 

END  IF 

IF  (J4.  EQ.  K)  THEN 
TEMP=I.  +TEMP 
END  IF 

H( J 4 , K ) =TEMP 
CONTINUE 
CONTINUE 

DO  131  J4=COLUM.N4 
DO  130  K=COLUM,N4 
TEMP=0. 

DO  140  L=COLUM.N4 

TEMP=TEMP+H(  J4, L) *A(  L, K) 

CONTINUE 
Cf J4,K)=TEMP 
CONTINUE 
CONTINUE 

DO  150  J4=C0LUM,N4 
TEMP=0. 

DO  160  L=COLUM(N4 

TEMP=TEMP+H(  J4 , L ) *F(  L ) 

CONTINUE 

T(J4)=TEMP 

CONTINUE 

DO  170  J4=COLUM,N4 
DO  180  K=COLUA,N4 

co«K,<,J4'K) 

F(  J4J=T(  J4) 

CONTINUE 

CONTINUE 

X(  N4 )  =F(  N4  ) /A(  N4,N4) 

I=N4- 1 
SUM2=0. 

M2=I+1 

DO  30  J4=M2.N4 

SUM2=SUM2+ A(  I , J4 ) *X(  J4 ) 

CONTINUE 


IF(  I.  GE.  1)  GO  TO  812 
DO  141  1=1. N1 
S2(  I)=X(l) 

CONTINUE 
DO  142  I =N3 , N4 
N5=I-N1 
Sl(  N5)=X( I) 

CONTINUE 
UP(  J)=0.  0 
DO  5001  1=1 ,  N 

UP( J) =UP^  J J  +  £ S2 <L7  j -C3  (  I ) ) *U(  J- 1 )  +  Sl(  I ) * Y(  J- 1 ) 
CONTINUE 

Uir5=jiiP(  J)+V<  J) )/( 1-S2(  1} } 

Y( J5 )=0. 0 
DO  30O3  1=1,  N 
M6=N1- I 

YrJ5T=Y(  J5)-C6(  I)*Y( J5-M6)+C5( I )*U(  J5-M6) 
CONTINUE 
DO  189  1=1 ,N1 
YM( I )=0. 0 


CONTINUE 
YM(  J5  )=0.  0 
DO  2568  1=1, N 
M8=N- I + 1 
Bl(  I  )=C5(  M8 ) 

CONTINUE 
DO  3004  1=1, N 

CONTI fjUE)=Y^(  J5)“C4<  I>*YM<  J5-I  )+Bl(  I )  *V(  J5-I) 

MK(  J5  )=-Y(  J5 ) 

DO  3065  1=1,  N 

CONTI^IUE }  =Mk(  J5  )+B{1)  *V<  J5‘ 1  >  "C4(  1  >  *Y(  J5‘ 1  > 
MU( J5)=ABS(  MK(J5 ) ) 

IF(MU(  J5).  LT.  20.  6)  THEN 

ELsic  )=^C(  J5) 

V(  J5)=VC(  J5 )+VP{  J5) 

END  IF 

DO  120  1=1, N2 
Fl(  I,M)=0.  0 
CONTINUE 
DO  118  1=1, N 
L1=N+ I 


J2=J+ 1- I 

Fl(  I,M)=-Y{  J2 ) 

FI  Ll,M)=U( J2 ) 

CONTINUE 
DO  119  1=1, N2 
FT(M,  I)=0.  0 
CONTINUE 

CALL  TRANS(  N2 ,M, I . K , FI , FT ) 

CALL  FAMILY? N2,M, *U17) 

CALL  RESULT? N2 , M , I , K , U1 ) 

DO  3001  I=1,N2 
U4( I,M)=6.  0 
CONTINUE 

CALL  SUB( N2 , N2 , I , L , U3 , U1 , U4 ) 

L6=J-N 

C7(  L6 ) =0. 0 

DO  3002  1=1 ,N2 

coN??^»-C’(L6)iU4(I-M)*U4(I-M) 

gSffi8RT(C7(L6,) 

WR I TE(  8 . 2  2  5 5 ) *  PARAMETER  ERROR 1 
FORMATS 12X.A15) 

DO  8882  1=1.100 

WRITER 8,8883)  C8(I) 

FORMAT( 15X,F12. 6) 

CONTINUE 


WRITE(8. 1111) ' SYS.  OUTPUT' , ' 
FORMAT (  9X . A1 1 ,6X , A1 1 ) 

DO  976  I=1.106 

WRITE  (6, 8602 )  Y(  I ) , YM(  I ) 

WRITE  (8.8002)  Y?I),YM(I 

FORMAT?  9&,  F10.  6 , 6X,  fclO;  6 ) 

CONTINUE 


MOD.  OUTPUT' 


WRITE  (8.8002)  Y(I),YM(I) 
FORMAT?  9&,  F10.  6 , 6X,  fclO;  6 ) 
CONTINUE 
DO  1025  1=1, L8 
SE( I,l)«fa I) 

SE? I.2)=YM(  I) 

KD  l)=I 
CONTINUE 
C8(  36 )=0. 65827 
C81  37  =0.  65827 
C81  38  =0.  65827 
C8  39  =0. 65827 
C81  40  =0.  65827 
CALL  VPLOT(KD, SE.L8, 2 ,  TIME 
CALL  V1FL0T( KD, C8, L8 , 1, ’TIME 
CALL  V2PL0Ti  KD, MU,  L8, 1 .  TIMI 
CALL  V3PL0Ti  KD,V,L8,i,  TIME 
CALL  V4PL0TI1  KD,U,  L8, 1,  TIME 


VPLOTfKD, SE.L8, 2 ,  TIME  '  ' , 
VIPLOTf  KD,C8, L8, 1, ;TIME  *  ' 
V2PL0T?  KD, MU, L8, 1 .  TIME  '.' 


OO  t-KSCO  if)  & 

Ot^CO  O  OXT'OI  U  UCn  OO'  oo>  OO 


?  —  1  i  H 

8, 70)(  T(  P  .  S) ,  S=1  .  °) 

7T2;ii'  io(15x<  fi2. 


END 

♦THIS  SUBROUTINE  PRINTS  THE  ARRAYS  AS  AN  MATRIX  FORM. 
SUBROUTINE  RESULT( H, O, P , S, T) 

INTEGER  H  ,  0  ,  P  ,  S 
REAL  T  (  H  ,  0  ) 

DO  80  P  =  1  ,  H 
WRITE( 

PRINT 

F0RMAT(  T6'  7^2^  9%.',  10fSx*F12.  6) ) 

CONTINUE 
RETURN 
END 

***THIS  SUBROUTINE  CALCULATES  THE  MULTIPLICATION  OF 
TWO  MATRICES.  *** 

SUBROUT  I NE  C  ALCUL  l\JjV,Y,Z,XtZXtZX,UeQ) 

INTEGER  U ,V .YjZ.XjZX 
REAL  ZT(U/vK*(6#'fc)#Q(Y/V) 

DO  93  Z=1,U 
DO  92  X=1,V 
ZT(Z.X)  =0.0 
DO  91ZX=1,Y 

C0^)=6t(Z'X)+W(  Z/ZX)*c(  zx,x) 

CONTINUE 

CONTINUE 

RETURN 

END 

***THIS  SUBROUTINE  WRITES  THE  MATRIX  NAME  AND  ROW, 
COLUMN  NUMBERS.  ***/ 

SUBROUT I NE  FAM I LY(  ROW , COLUMN  #  NAME ) 

INTEGER  ROW, COLUMN 
CHARACTER* 2  NAME 
WRITE( 8, 94)  ' MATRIX  , NAME 
PRINT  94. 'MATRIX' .NAAe 

FORMATfTi , ' 0 ; , T2 , iOX, T12 . A6 , T19 , IX, T20 , A2 ) 

WRITE( 8,95)  ROW  NUMBER.  1 , ROW 
PRINT  96/ROW  NUMBER.  T. ROW 

FORMATS Ti , 1  0  ,  T2 . IOX , Ti 2 ,A1 1 , T23 . 8X , T3 1 , I 2 ) 

WRITE( o,96)  COLUMN  NUMBER.  T, COLUMN 

PRINT  96. rCOLUMN  NUMBER.  T, COLUMN 

FORMAT( Tl, r0  , T2 , IOX, T12 , A14, T26, 5X, T31 , 12) 

RETURN 

END 

***THIS  SUBROUTINE  CALCULATES  THE  SUMMATION  OF  THE 
absolute  values  OF  THE  SAME  COLUMN  ELEMENTS. *** 
SUBROUTINE  ALI{ II, 12 ,H1, H2 , Gl, W2 ) 

INTEGER  II, 12, HI. 62 
REAL 


DO  42  Hl=l , I 1 

„rtVtlT,TXTtTT,  W2(  1 ,  H2  )=&2(  1 ,  H2  )  +ABS(  Gl(  HI ,  H2  ) ) 

42  CONTINUE 

41  CONTINUE 

RETURN 

END 

C  ***THIS  SUBROUTINE  MULTIPLIES  THE  MATRIX  BY  SCALAR 
NUMBER.  *** 

SUBROUTINE  SCALAR( K1 , K2 , Gil , G2 , PI , G3 , G4) 

INTEGER  K1,K2,G11.G2 
REAL  G3(Kl,K2j  . G4( K1 , K2 ) , PI 
DO  43  G2=1,k2 

DO  44  G11=1,K1 

G4( G11,G2)=G3(  G11,G2)*P1 

44  CONTINUE 

43  CONTINUE 
RETURN 
END 

C  ***THIS  SUBROUTINE  CALCULATES  THE  SUM  OF  THE  TWO 
MATRICES.  *** 

SUBROUTINE  SUM( K3.K4, 15 , K5 , XI , X2 , X3 ) 

INTEGER  K3,K4,I5.fc5 

REAL  X1^K3  ^|4^X2  (  K3  ,  K4 ) ,  X3  (  K3 ,  K4 ) 

DO  46' I 5=1. K3 

X3( I5,K5)=X1( I5,K5)+X2(  15, K5) 

46  CONTINUE 

45  CONTINUE 
RETURN 
END 

C  ***THIS  SUBROUTINE  COPPIES  THE  MATRICES.  *** 

SUBROUT I NE  COP YM(  K3 , K4 , 1 5 , K5 , XI , X2 ) 

INTEGER  K3<K4.I5'K5 
REAL  Xl(K3,K4j ,X2(K3,K4) 

DO  45  X$=1,K4 

DO  46  15=1, K3 
Xl(  I5,K5)=fc2{  1 5 ,  K5 ) 

46  .  CONTINUE 

45  CONTINUE 
RETURN 
END 

C  *THIS  SUBROUTINE  FINDS  THE  TRANSPOSE  OF  THE  MATRIX.* 
SUBROUTINE  TRANS(  K3 , K4 , 1 5 , K5 , XI , X2 ) 

INTEGER  K3,K4.I5.K5 
REAL  X1(K3 , K4) ,  X2(  K4, K3 ) 

DO  45  k£=1,K4 

DO  46  1 5=1 , K3 
X2(  K5 ,  I5)=&1(  15, K5 ) 

46  CONTINUE 

45  CONTINUE 


CONTINUE 

CONTINUE 


RETURN 

END 

C  ***THIS  SUBROUTINE  CALCULATES  THE  SUB.  OF  THE 
TWO  MATRICES.*** 

SUBROUTINE  SUB(K3.K4, 15 , K5 , XI , X2 , X3 ) 

INTEGER  K3 , K4, 15 ,  K5 

REAL  X1(K3 , K4  j ,  X2  ( K3 , K4 ) , X3  (  K3 , K4 ) 

DO  45  K5 — I , K4 

DO  46  15=1. K3 

X3( I5,K5)=X1( I5/K5)-X2(  15, K5) 
46  CONTINUE  ' 

45  CONTINUE 

RETURN 
END 

0  ****************************************** 

C  *  THIS  SUB.  PLOTS  THE  ARRAYS  USING  * 

C  *  THE  DISSPLA  PROGRAM.  * 

0  ****************************************** 

C 

SUBROUTINE  VPLOT(  T, U, M, N, LABELX, LABELY) 
0*****.********  '  ' 
REAL  U(800, 4)  ,  T(800 ) ,  f(  800) , YX(  4)  ,  YN(  4) 
INTEGER  PMAX, PMIN, M,N,  LABEwC,  LABELY 


CALL  AMAX(  T , M , TMAX , PMAX ) 
CALL  AMIN(T, M,TMIN, PMIN) 
DO  20  J  =1,6 
DO  10  I  =1 , M 


CONT 


He=",wi 


CALL  AMAXf Y. M, YMAX, PMAX) 

YY/  Tl  —  VM&Y 

CALL 'aMIN(Y,M,YMIN, PMIN) 
YN( J)  =YMIN 
WRITE( 6,11)  YMAX, YMIN 
CONTINUE 


CALL  AMIN 


WRITE( 6, 11) 

CALL  TEfc  6 18 
CALL  BLOWUP  {  1.2) 

CALL  PAGE( 11.  ,8.  5) 
CALL  NOBRDR 
CALL  AREA2D(  9.  ,6.  ) 
CALL  XNAME( LABELX, 6) 
CALL  YNAME(  LABELY,  6) 
CALL  HEAD IN  ( tPLANT 
CALL  CROSS 


X{  YX,N,YP1 
N(  YN,  N,  YW 

ll)  ItmAx, 


YMIN, PMIN 
X,YMIN 


AND  MODEL  OUTPUTS  $',100,1. ,1) 


CALL  GRAF( TMIN,  SCALE  , TMAX, YMIN, 'SCALE' , YMAX) 
CALL  SPLINE 
DO  40  J  =1 , N 
DO  30  I  =1 . M 
Y[lj  =U(l,J) 

30  CONTINUE 

IF(J.  EQ.  2)  CALL  CHNDOT 
CALL  CURVE  (  T,  Y,  M,  0 ) 

40  CONTINUE 

CALL  ENDPL(O) 

C  CALL  DONEPL 

11  FORMAT(  ,  1 , 10X,2G12.  5/) 

13  FORMAT  '  ' , 10X,4G12.  5/) 

RETURN 

END 

C 

SUBROUTINE  V1PL0T(  T,U,M,N, LABELX, LABELY) 

C  REAL  U(800,4) ,T( 800) ,Y( 800) ,YXf 4).YN(  4) 

INTEGER  PMAX, PMIN, M,N, LABELX,  IABELY 
C 

CALL  AMAXf  T , M , TMAX , PMAX ) 


WAX(T,M, 

call  amin[t,m,tmin;pmin) 

DO  20  J  =l,ft 
DO  10  I  =1,M 
Y(  I )  =U(I,J) 

CONTINUE 

CALL  AMAXJY.M, YMAX, PMAX) 

VY l  =  VM1Y 

CALL  AMIN(Y,M, YMIN, PMIN) 


TMIN, PMIN] 


CALL 
YN(  j: 


uyiiN(  y, 

=YMIN 


WRlTE(  6,11)  YMAX, YMIN 
CONTINUE 


X(YX,N, YMAX, PMAX) 

CALL  AMIN(YN,N. YMIN, PMIN) 

WRITE(6. 11)  £mAx,YMIN 
CALL  TEfc  6l8 
CALL  BLOWUP( 1.2) 

CALL  PAGE(  11.  ,8.  5) 

CALL  NOBRDR 
CALL  AREA2D(  9.  ,6.  ) 

CALL  XNAME(  LABELX, 6 ) 

CALL  YNAME( LABELY. 61 

CALL  HEAD IN  ( tPARAMETER  ERROR  $ ' , 100, 1. , 1 ) 

CALL  CROSS 

CALL  GRAF( TMIN, 'SCALE' , TMAX, YMIN, 'SCALE' , YMAX) 
CALL  SPLINE 
DO  40  J  =1 , N 


AREA2D(  9.  ,6.  ) 
XNAME( LABELX, 6) 
YNAME  LABELY. b) 
HEADIN  (tPARAME 


TER  ERROR  $',100,1. ,1) 


DO  30  I  =1,M 


30 

CONTI: 

IF(  J. 

CALL  i 

40 

CONTINUE 

CHNDOT 

M,0) 


CALL  ENDPL(O) 

CALL  DONEPL 

FORMAT (  ,  ;/10X,2G12. 

FORMAT  '  , 10X, 4G12. 5 

RETURN 

END 


SUBROUTINE  V2PL0T( T,  U,M, N, LABELX, LABELY) 

C  REAL  U( 800,4) ,T(800) , Y( 800)  YXf 4) ,YN( 4) 
INTEGER  PMAX , PMI N ,  M ,  N , LABELX , LABELY 
C 

CALL  AMAX( T ,  M , TMAX , PMAX ) 

CALL  AMINi T, M, TMIN, PMIN) 

DO  20  J  =1,N 
DO  10  I  =1 . M 
Y(  I )  =U(I,J) 

10  CONTINUE 

CALL  AMAXf Y . M y YMAX , PMAX ) 

VY /  T\  —  VMay 

CALL  AMIN(  Y, M, YMIN, PMIN) 

YN(  J)  =YMIN 
WRITE(  6,11)  YMAX, YMIN 
20  CONTINUE 

WRITEf  6, 13 )(  YX(  I ) , 1=1 , 4) 


YMIN, PMIN 
X,YMtN 


WRITE(  6 . 13  H  YN(  Ij  ,  1=1, 4j 
CALL  AmAx(  YX,N, YMAx,PMAX) 

CALL  AMIN(YN,N. YMIN, PMIN) 

WRITE(6, 11)  -$rMAx,YMiN 
CALL  TE k  618 
CALL  BLOWUP( 1.2) 

CALL  PAGE( 11.  ,8.  5) 

CALL  NOBRDR 
CALL  AREA2D(  9.  ,6.  ) 

CALL  XNAME(  LABELX, 6) 

CALL  YNAME  LABELY, 61 

CALL  HEADIN  ( T OUTPUT  PREDICTION  ERROR  $’,100,1. 
CALL  CROSS 

CALL  GRAF( TMIN, ’SCALE’ , TMAX, YMIN, ’SCALE’ , YMAX) 
CALL  SPLINE 
DO  40  J  =1 , N 
DO  30  I  =1 , M 


rtriEirtZUI  -a.  .  o.  i 

XNAME(  LABELX, 6) 
YNAME  LABELY, 61 
HEADIN  (TOUTPUT 


PREDICTION  ERROR  $’,100,1. ,1) 


IF(J.  EQ.  2)  CALL  CHNDOT 
CALL  CURVE  (T,Y,M,0) 
CONTINUE 
CALL  ENDPL(O) 

CALL  DONEPL 

FORMAT (  ,  ; , 10X,2G12.  5/) 
FORMAT?  ' ,10X,4G12.  5/) 
RETURN 


SUBROUTINE  V3 PLOT(  T , U , M , N , LABELX , LABELY ) 
0** *********** 

REAL  U(800, 4) ,  T(800 ) ,  Y( 800 ) , YX( 4) , YN( 4) 
INTEGER  PMAx,  £MIN,M,  fi,LABElJc,  LABELY 


CALL  AMAXf  T ,  M , TMAX , PMAX ) 

CALL  AMIN? T, M,TMIN, PMIN) 

DO  20  J  =l,ft- 
DO  10  I  =1 , M 

CONT-tN1jrE-U  I/J) 

CALL  AMAJKY.M,  YMAX,  PMAX) 

YX( J)  =  YMAX 

CALL  AMIN(Y,M, YMIN, PMIN) 

YN(J)  =YMIN 
WRITE( 6,11)  YMAX, YMIN 
CONTINUE 

WRITE(6,13)(YX(  I),  1=1, 4) 

WRITE?  6, 13  j?  YN?  IJ  .  1=1 . 4 j 
CALL  AMAX(YX,N, YMAX, PMAX) 

CALL  AMIN? YN,N. YMIN, PMIN 
WRITE(6, 11 J  YMAX, YMIN 
CALL  TEK  618 
CALL  BLOWUP( 1.2) 

CALL  PAGE?*  11.  ,8.  5) 

CALL  NOBRDR 
CALL  AREA2D( 9.  , 6.  ) 

CALL  XNAMEf  LABfeLX ,  6 ) 

CALL  YNAME?  LABELY, 6) 

CALL  HEADIN  ( tEXTERNAL  INPUT  $',100,1. ,1) 

CALL  CROSS 

CALL  GRAF(TMIN, 'SCALE' , TMAX, YMIN, 'SCALE' , YMAX) 
CALL  SPLINE 
DO  40  J  =1 , N 
DO  30  I  =1 . M 
Y(IJ  =U?I,J) 

CONTINUE 

IF(J.  EO.  2)  CALL  CHNDOT 
CALL  CURVE  (T,Y,M,0) 

CONTINUE 


TMAX, PMAX) 
TMIN, PMIN) 


AREA2D(  9.  ,6.  ) 
XNAMEf LAB&LX, 6) 
YNAME  LABELY, 6J 
HEADIN  (tEXTERN, 


IAL  INPUT  $'  ,100,1.  ,1) 


V'-.V-w* 


•;o 


CALL  ENDPL(O) 

C  CALL  DONEPL 

11  FORMAT (  ^  : ,10X,2G12.  5/) 

13  FORMAT ( '  ' , 10X,4G12.  5/ 

RETURN 

END 

SUBROUTINE  V4PLOT(  T ,  U ,  M ,  N , LABELX , LABELY ) 
q************* 


REAL  U(  800 ,  4 ) ,  T{  800 ) , Y(  800 ) . YX(  4 ) , YN( 4 ) 
INTEGER  PMAX, PMIN, M, N, LaBEla, LABELY 

CALL  AMAX(  T,M, TMAX, PMAX) 

CALL  AMIN( T, M, TMIN, PMIN) 

DO  20  J  =1,6 


CALL  AMAX(  T,M, TMAX, PMAX) 
CALL  AMIN( T, M, TMIN, PMIN) 

DO  20  J  =1,6 
DO  10  I  =1,M 

cont^n^e-^  1 ' 

CALL  AMAXf Y.M, YMAX, PMAX) 
YX( J)  =  YMA& 

CALL  AMIN(Y,M,YMIN, PMIN) 
YN( J)  =YMIN 
WRITE( 6, 11 )  YMAX, YMIN 
CONTINUE 

WRITE(  6, 13 )  (  YX(  I), 1=1, 4) 
WRITE(  6 . 13  j  (  YN(  I) ,1=1, 4 j 
CALL  AmAX(YX,N, YMAX, PMAX) 
CALL  AMINfYN.N. YMIN, PMIN) 
WRITE (6, ll )  YMAX, YMIN 
CALL  TE6  618 
CALL  BLOWUPf  1.2) 

CALL  PAGEfll.  ,8.  5) 


CALL  NOBRDR 
CALL  AREA2D(  9.  ,6.  ) 

CALL  XNAME(LAB6LX,6) 

CALL  YNAMEl  LABELY, 6 ) 

CALL  HEADIN  riNP6T  TO  THE  PLANT  $',100,  1.  ,1) 
CALL  CROSS  .  , 

CALL  GRAF(  TMIN, ’ SCALE ' , TMAX, YMIN, ' SCALE ' , YMAX) 
CALL  SPLINE 
DO  40  J  =1 , N 
DO  30  I  =1 - M 

IFf  J.  EQ.  2)  CALL  CHNDOT 
CALL  CURVE  (T,Y,M,0) 

CONTINUE 
CALL  ENDPL(O) 

CALL  DONEPt 

FORMATS ,  , , 10X, 2G12.  5/) 

FORMAT(  ’  ' ,10X,4G12.  57) 


FORMATi 


,  10X, 4G12. 5 


v  »  *  w~„  r .  r  wv  *~  «~  v'  *r*  »*  „■ 
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